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Abstract 

This  research  effort  seeks  to  characterize  a  vision-aided  approach  for  an  Un¬ 
manned  Aerial  System  (UAS)  to  autonomously  determine  relative  position  to  another 
aircraft  in  a  formation,  specifically  to  address  the  autonomous  aerial  refueling  prob¬ 
lem.  A  system  consisting  of  a  monocular  digital  camera  coupled  with  inertial  sensors 
onboard  the  UAS  is  analyzed  for  feasibility  of  using  this  vision-aided  approach. 

A  three-dimensional  rendering  of  the  tanker  aircraft  is  used  to  generate  predicted 
images  of  the  tanker  as  seen  by  the  receiver  aircraft.  A  rigorous  error  model  is 
developed  to  model  the  relative  dynamics  between  an  INS-equipped  receiver  and  the 
tanker  aircraft.  A  thorough  image  processing  analysis  is  performed  to  determine  error 
observability  between  the  predicted  and  true  images  using  sum-squared  difference  and 
gradient  techniques.  To  quantify  the  errors  between  the  predicted  and  true  images, 
an  image  update  function  is  developed  using  perturbation  techniques.  Based  on  this 
residual  measurement  and  the  inertial/dynamics  propagation,  an  Extended  Kalman 
Filter  (EKF)  is  used  to  predict  the  relative  position  and  orientation  of  the  tanker 
from  the  receiver  aircraft. 

The  EKF  is  simulated  through  various  formation  positions  during  typical  aerial 
refueling  operations.  Various  grades  of  inertial  sensors  are  simulated  during  different 
trajectories  to  analyze  the  system’s  ability  to  accurately  and  robustly  track  the  relative 
position  between  the  two  aircraft. 
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Using  Predictive  Rendering  as  a  Vision-Aided  Technique 
for  Autonomous  Aerial  Refueling 


I.  Introduction 

The  advent  of  unmanned  aerial  systems  (UAS)  has  revolutionized  the  way  the 
LInited  States  Air  Force  (USAF)  conducts  operations.  No  longer  are  aerial 
platforms  limited  by  human  endurance,  but  rather  the  endurance  of  the  aircraft  itself. 
This  limiting  factor  can  be  even  further  reduced  with  the  use  of  automated  in-flight 
LIAS  refueling.  The  first  unmanned  systems  (RQ-4  Global  Hawk  and  MQ-1  Predator) 
were  used  exclusively  for  reconnaissance  and  could  loiter  over  a  target  for  extreme 
lengths  of  time.  The  Predator  soon  transitioned  into  an  armed  reconnaissance  role. 
Today’s  systems  have  evolved  even  further  into  hunter/killer  systems  such  as  the  MQ- 
4  Reaper.  One  can  easily  see  this  trend  will  develop  unmanned  systems  into  a  more 
complex,  larger  component  of  the  battlefield  in  the  years  to  come  [12], 

With  the  imminent  arrival  of  Unmanned  Combat  Aerial  Systems  (UCAS),  the 
need  for  tactical  refueling  will  be  vital  for  mission  success.  Designed  for  high  perfor¬ 
mance  and  larger,  more  complicated  weapon  payloads,  these  systems  will  have  both 
shorter  range  and  loiter  time  than  the  current  unmanned  systems  fielded  today.  These 
systems  will  likely  be  a  critical  component  of  the  tactical  battlefield,  and  comman¬ 
ders  will  require  them  to  perform  missions  that  are  currently  dominated  by  manned 
systems.  Thus,  the  need  for  a  robust,  autonomous  refueling  technique  is  vital  for 
force  projection  of  unmanned  systems  and  overall  mission  success  in  future  USAF 
operations  [2,8,12], 

Many  efforts  have  already  shown  success  in  solving  this  problem  (articulated  in 
detail  in  Chapter  II),  primarily  using  Global  Positioning  System  (GPS)  approaches  [8, 
13].  With  GPS,  possible  satellite  acquisition  issues  exist  for  the  receiver  aircraft  due 
to  the  tanker  blocking  large  areas  of  the  sky  above  it  [9].  Furthermore,  relying  on 
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GPS  solutions  in  every  combat  environment  could  be  detrimental  when  GPS  denial 
and  degradation  exist  due  to  jamming  or  sensor  failure.  The  Air  Force  Research 
Laboratory’s  (AFRL)  Automated  Aerial  Refueling  (AAR)  Program  has  an  interest 
in  determining  how  well  an  autonomous  system  can  function  in  environments  where 
GPS  is  either  degraded  or  denied  [20].  This  research  will  directly  help  analyze  this 
specific  area  of  the  AAR  program. 

1.1  Methods  of  Aerial  Refueling 

There  are  two  standard  methods  for  in-flight  refueling  operations.  The  first, 
employed  by  the  United  States  Navy  and  Marine  Corps,  is  known  as  the  probe  and 
drogue  technique.  For  this  method,  the  refueling  aircraft  has  a  refueling  probe  that 
makes  contact  with  a  drogue  chute  on  a  refueling  hose  that  trails  from  the  tanker 
aircraft.  Figure  1.1  shows  a  Marine  Corps  F/A-18  Hornet  receiving  fuel  from  an  Air 
Force  KC-10  Extender. 

The  second  major  method  of  aerial  refueling  is  the  flying  boom  method  used  by 
the  USAF.  For  this  approach,  the  receiving  aircraft  maintains  a  formation  position 
just  below  the  tail  of  the  tanker  aircraft,  and  the  boom  operator  onboard  the  tanker 
inserts  the  refueling  boom  into  the  refueling  receptacle  of  the  receiver.  Figure  1.2 
displays  this  technique. 


Figure  1.1:  An  F/A-18  Hornet  receiving  fuel  from  a  KC-10 
Extender  using  the  probe  and  drogue  aerial  refueling  method 
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Figure  1.2:  An  F-22  Raptor  receiving  fuel  from  a  KC-135 

Stratotanker  using  the  aerial  boom  refueling  method 

1.2  Research  Objectives 

By  applying  an  image  processing  algorithm  to  images  received  from  a  monocular 
camera,  this  research  will  investigate  if  the  receiver  aircraft  can  accurately  determine 
the  relative  attitude  and  range  to  the  lead  tanker  aircraft.  Fusing  this  image  data  with 
the  onboard  Inertial  Navigation  System  (INS)  data  will  provide  a  robust  navigation 
technique  for  autonomous  UAS  formation  flight  that  could  potentially  be  independent 
of  external  sources  such  as  GPS.  The  ultimate  goal  will  be  validating  a  robust,  vision- 
aided  technique  for  accurate  relative  positioning  between  the  tanker  and  receiver 
aircraft. 

1 . 3  Scope 

This  research  will  focus  on  achieving  accurate  relative  position  of  a  receiver 
aircraft  to  a  KC-135R  Stratotanker.  As  aforementioned,  the  baseline  system  will 
have  a  monocular  camera  as  well  as  data  from  an  on-board  INS.  The  receiver  aircraft 
is  assumed  to  start  in  a  position  behind  and  below  the  tanker,  using  the  boom  method 
of  aerial  refuclng.  This  will  ensure  the  images  are  free  from  unwanted  objects  such 
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as  ground  clutter  and  other  aircraft.  This  research  will  analyze  potential  differences 
in  INS  quality  and  consider  various  image-processing  techniques  to  develop  a  robust, 
accurate  approach. 

1.4  Overview 

This  thesis  will  be  organized  as  follows.  Chapter  II  will  highlight  key  math¬ 
ematical  background  to  include  applicable  coordinate  systems  used  in  the  relative 
navigation  algorithm  as  well  as  transformation  techniques  between  them.  Kalman 
filtering  will  be  discussed  in  detail,  to  include  the  basic  Kalman  Filter  as  well  as  the 
specific  implementation  used  in  this  research,  the  Extended  Kalman  Filter  (EKF). 
Then,  fundamentals  of  an  INS  will  be  discussed  as  well  as  example  mechanization 
equations.  Also,  some  fundamental  concepts  of  imaging  sensors  are  presented.  The 
chapter  concludes  by  analyzing  previous  work  concerning  the  autonomous  aerial  refu¬ 
eling  problem,  and  how  this  research  seeks  to  complement  those  previous  efforts  and 
exploit  any  lessons  learned. 

Chapter  III  presents  the  development  of  a  predictive  rendering  environment  to 
generate  predicted  images  of  the  tanker  aircraft  as  seen  by  the  receiver.  Several  image- 
processing  techniques  (sum-squared  difference,  magnitude  of  gradient,  and  magnitude 
of  gradient  with  threshold)  are  analyzed  to  determine  their  suitability  for  detecting 
errors  between  the  true  and  predicted  images.  A  relative  dynamics  model  is  developed 
to  simulate  the  wing  and  lead  aircraft  flight  characteristics  during  the  refueling  oper¬ 
ation.  A  discussion  of  the  EKF  implementation,  including  the  use  of  image  updates, 
concludes  the  chapter. 

Chapter  IV  analyzes  the  results  of  various  simulations  using  the  filter  described 
in  Chapter  III.  Several  trajectories  are  analyzed  with  the  receiver  and  tanker  aircraft 
in  different  positions  of  typical  refueling  formations  (e.g.,  pre-contact  and  contact). 
A  presentation  of  filter  performance  with  various  grades  of  inertial  sensors  and  initial 
errors  is  also  presented  based  on  the  simulation  results. 
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Finally,  Chapter  V  will  highlight  the  overall  conclusions  of  the  image-aided 
filter  for  autonomous  refueling,  analyzing  the  overall  performance  of  the  technique. 
An  overall  analysis  of  the  research  accomplishments  and  recommendations  for  future 
research  conclude  the  chapter. 
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II.  Mathematical  Background  and  Previous  Research 


This  chapter  describes  the  key  mathematical  and  technical  background  infor¬ 
mation  necessary  to  understand  the  algorithm  development  of  Chapter  III. 
These  key  concepts  are  highlighted  in  the  beginning  sections  of  this  chapter  and 
include:  Kalman  filtering,  coordinate  systems  and  transformations,  INS  description 
and  mechanization,  optical  sensors,  and  image  processing.  The  last  section  of  the 
chapter  highlights  the  current  state  of  autonomous  aerial  refueling  research. 

2.1  General  Mathematical  Notation 

Before  developing  mathematical  equations,  a  brief  description  of  the  general  no¬ 
tation  is  in  order.  For  this  document,  scalar  quantities  will  be  identified  by  lower-case 
variables  (e.g.,  x ).  Vector  quantities  will  be  identified  by  bold,  lower-case  variables 
(e.g.,  x),  and  matrices  will  be  identified  by  bold,  upper-case  variables  (e.g.,  X). 

2. 2  Kalman  Filtering 

The  Kalman  Filter  is  an  optimal  estimation  technique  frequently  employed  in 
navigation  systems  that  uses  statistical  techniques  to  estimate  the  navigation  state. 
There  are  many  varieties  of  this  technique,  but  the  two  discussed  here  are  the  tra¬ 
ditional  Kalman  Filter  and  the  Extended  Kalman  Filter.  Before  constructing  this 
optimal  estimator,  a  mathematical  model  must  be  constructed  to  approximate  the 
dynamics  of  the  system.  Using  a  state-space  format,  the  state  vector  of  interest,  x, 
changes  over  time  by  the  following  expression 

x(f)  =  F(t)x(f)  +  B(f)u(f)  +  G(f)w(f)  (2.1) 

where  F(f)  is  the  state  dynamics  matrix,  B(f)  is  the  input  matrix,  u(f)  is  the  vector 
of  deterministic  inputs,  G  (f)  is  the  noise  distribution  matrix,  and  w (f)  is  the  noise 
vector.  The  values  in  the  noise  vector  are  typically  assumed  to  be  white,  Gaussian 
noise  (WGN)  sources.  Similarly,  discrete-time  measurements  are  used  to  update  the 
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filter  and  are  typically  of  the  form: 


z  (tfc)  =  H(ifc)x(tfc)  +  v(4)  (2.2) 

where  z(4)  are  the  discrete  time  measurements,  H(4)  is  the  measurement  observ¬ 
ability  matrix,  and  v(4)  is  the  discretized  measurement  noise  vector. 

2.2.1  Traditional  Kalman  Filter.  Now,  with  a  stochastic  model  of  the  sys¬ 
tem,  the  traditional  Kalman  Filter  can  be  constructed.  This  estimator  tracks  a  state 
estimate  (x)  and  a  state  error  covariance  (P)  through  an  alternating  series  of  time- 
propagation  and  measurement  updates.  These  and  other  variables  used  in  the  filter 
will  often  have  a  subscript  and  superscript  notation.  The  subscript  notation  indicates 
at  which  time  epoch  the  variable  is  evaluated.  The  superscript  indicates  whether 
this  is  immediately  before  or  after  a  measurement  update  occurs.  For  example,  the 
notation  of  PjT  indicates  the  a  priori  state  error  covariance  estimate  just  before  a 
measurement  is  taken  at  some  time  epoch  k.  The  notation  of  P^  indicates  the  a 
posteriori  state  error  covariance  estimate  just  after  a  measurement  is  taken  at  some 
time  epoch  k. 

To  construct  the  filter,  the  first  step  is  to  write  the  discrete-time  stochastic 
difference  equation  from  the  continuous  time  model  shown  above.  Before  writing  this 
equation  directly,  there  are  a  few  specific  calculations  that  must  be  made.  The  needed 
values  are  B^,  $(A t),  and  Q/,..  is  the  discretized  version  of  the  input  matrix  B(t), 
$(A t)  is  the  discrete  time  state  transition  matrix  over  the  specific  time  interval,  and 
Qfc  is  the  discretized  process  noise.  These  matrices  are  calculated  as  follows  [11]: 

rtk+i 

Bk  =  /  $(4+i  -  r)Bu(r)dr  (2.3) 

Jtk 

where: 

$(A  t)  =  eFAi  (2.4) 
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Qk  is  solved  using  the  Van  Loan  calculation,  which  is  performed  as  follows: 


A*  = 


-F  GQG1 
0  FT 


B*  =  eA*  = 

1 

a 

H0H 

Bn 

l - 

CM 

m 

0  <FT 

Cd 

to 

B22  I 

Q/t  —  B22B12 

Now,  the  stochastic  difference  equation  can  be  formed  as  seen  in  Equation  2.8: 


(2.5) 

(2.6) 

(2.7) 


X;.  I  -  +  11;  11/  +  w; 


(2.8) 


Now  ail  the  necessary  elements  exist  to  write  the  equations  for  executing  the  Kalman 
filter.  The  first  function  of  the  Kalman  filter,  time  propagation,  is  accomplished  with 


the  following  two  calculations  for  each  time  step  [11]: 

Xfc  =  +  Bfc.iUfc.!  (2.9) 

Pfc  =  &k-iPk-i*k-i  +  Qfc-i  (2-10) 

Then,  at  each  time  epoch  k  where  a  measurement  is  available,  the  following  measure¬ 
ment  update  calculations  are  made: 

K*  =  P*  +  R*]-1  (2.11) 

=  +Kr(zt-Hrxj)  (2.12) 

Pt+  =  Pr  -  K,  H,P^  (2,13) 


The  matrix  K/,  is  known  as  the  Kalman  gain  and  the  matrix  R*,  is  the  measurement 
uncertainty.  This  calculation  process  is  repeated  for  the  duration  of  filter  execu¬ 
tion  [11], 


2.2.2  Extended  Kalman  Filter.  In  many  system  models,  the  system  dy¬ 
namics  and/or  measurements  are  nonlinear  functions,  and  cannot  be  simply  written 
in  terms  of  matrix  multiplications  as  seen  in  Equations  2.1  and  2.2.  The  nonlinear 
dynamics  and  measurements  are  typically  written  in  the  following  form  [10]: 

x(t)  =  f  [x(t),  u(t),t]  +  G(t)w(t)  (2.14) 

z (tj)  =  h  [x(U)]  +  v(ti)  (2.15) 

The  Extended  Kalman  Filter  (EKF)  is  one  mechanism  to  account  for  these  nonlinear 
effects.  The  EKF  is  an  iterative  algorithm  where  the  filter  is  relinearized  about  each 
state  estimate,  x(t/~),  once  it  is  computed.  This  results  in  creating  a  new,  updated 
state  trajectory  in  order  to  better  estimate  the  next  time  propagation.  This  process 
should  minimize  the  difference  between  the  nominal  and  true  trajectories,  and  hence 
improve  the  overall  model  [10]. 

In  the  general  case,  at  a  time  U,  immediately  after  a  measurement  is  incor¬ 
porated  and  a  relinearization  about  the  new  state  estimate  is  calculated,  the  best 
estimate  of  5x(U)  is  5x(t?+)  =  0.  This  is  due  to  the  relinearization  just  accomplished 
about  the  best  state  estimate  [10].  Since  the  current  estimate  is  zero,  propagating  it  to 
the  next  sample  time  ti+ 1  will  also  yield  zero.  Using  this  knowledge,  the  measurement 
update  for  the  rclinearized  filter  is: 


6-x.(tf+1)  =  K(U+i)(zj+i  -  h[xn(ti+i/ti),ti+1\)  (2.16) 

The  notation  x„(U+i/U)  is  used  to  define  the  nominal  state  over  the  interval  tt  to  tl+\. 

Now  to  find  the  optimal  estimate  of  the  full  state  x(f/t,;),  simply  combine  the 
nominal  with  the  perturbation  estimate  [10]: 

x(t/U)  =  x„(t/U)  +  6x(t/ti )  (2.17) 
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The  equations  that  follow  in  this  section  form  the  iterative  algorithm  for  executing 
the  EKF.  With  the  given  dynamics  model  in  the  previous  section,  the  measurement 
update  step  at  time  tt  incorporates  the  current  measurement,  z,,  through  the  follow- 
tag  [10]: 

K(ij)  =  P(«r)HT[ii;xW)](H[(j;x(«r)]P(t-)HT[ij;x(«r)]  +  R(ti))-‘  (2.18) 

*(tf)  =  *(*,“)  +  K(*i)(zi  -  (2.19) 

P(«+)  =  P(i-)  -  K(t,)H[ti;x(i-)]P(f-)  (2.20) 

where  H[fj;x(f“)]  is  defined  as  the  partial  derivative  matrix  [10]: 

(2.21) 

x=x(q-) 

The  estimate  is  propagated  to  the  next  sample  time  fi+1  by  integrating  the  following 
differential  equations: 

Mt/U)  =  mt/U),  u(t),  t]  (2.22) 

P( t/U)  =  F[t;  k(t/ti)]P(t/ti)  +  P(t/t,)FT[f,  k(t/U )]  +  G(t)Q(t)GT(t)  (2.23) 

The  initial  conditions  for  these  integrations  are  based  on  the  measurement  update 
values  in  Equations  2.19  and  2.20. 


H [U'Mu );  = 


d  /i[x,  U 


dx. 


Mti/t,)  =  ±(tf) 

(2.24) 

P(U/U)  =  P(«+) 

(2.25) 

The  matrix  F [t;x(t/tj)]  from  Equation  2.23  is  defined  as  [10]: 


F  [t;k(t/ti)'_  = 


<9/[x,u(f),t] 


<9x 


(2.26) 


x=ic(t/ti) 
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The  linear  stochastic  difference  equations  can  then  be  written  in  the  form: 


5x(ti+i)  =  &{ti+1,  ti)Sx(ti)  +  w  d(ti)  (2.27) 

P  (U)~  =  &{ti+ 1,  ti)P(ti)$T(ti+i,  U)  +  Q  d{ti)  (2.28) 

After  propagating  to  the  next  sample  time,  the  new  values  are  defined  as: 

x(«“+i)  =  x(«i+1/ii)  (2.29) 

P(f.-+i)  =  P(W<<)  (2.30) 


These  propagate  and  update  steps  are  iterated  throughout  the  execution  of  the 
EKF  [10]. 

2. 3  Reference  Frames 

Before  delving  into  the  details  of  inertial  navigation  systems,  one  must  first 
understand  the  reference  frame  about  which  navigation  systems  operate  and  define 
position.  This  research  will  primarily  utilize  frames  presented  in  [17].  All  are  orthog¬ 
onal  and  follow  the  right-hand  convention,  and  are  defined  in  the  paragraphs  that 
follow.  A  graphical  representation  from  [17]  is  shown  in  Figure  2.1. 

The  inertial  frame  ( I-frame )  has  fixed  axes  that  do  not  rotate  in  space  in  relation 
to  the  fixed  stars.  This  frame  is  where  Newton’s  laws  of  motion  are  valid. 

The  Earth-fixed  inertial  frame  ( i-frame )  has  origin  at  the  center  of  the  Earth 
with  the  z  axis  aligned  with  the  axis  of  rotation,  and  x  and  y  axes  in  the  equatorial 
plane.  The  axes  do  not  rotate  in  relation  to  the  fixed  stars. 

The  Earth-centered  Earth-fixed  frame  (e- frame)  also  has  origin  at  the  center  of 
the  Earth.  The  z  axis  is  aligned  with  the  axis  of  rotation,  but  the  x  is  aligned  with 
the  intersection  of  the  equatorial  plane  and  the  prime  meridian  (Greenwich).  The 
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Figure  2.1:  Diagram  of  Navigation  Reference  Frames  [17] 


y  axis  is  orthogonal  based  on  the  right-hand  convention.  The  e-frame  rotates  with 
respect  to  the  inertial  frame. 

The  navigation  frame  ( n-frame )  has  its  origin  at  the  location  of  the  inertial 
navigation  system  (INS)  or  a  fixed  point  on  the  surface  of  the  Earth.  The  x  and  y 
axes  are  parallel  to  the  surface  of  the  Earth  (aligned  with  North  and  East  respectively), 
and  the  z  axis  points  in  the  direction  of  local  gravity  (down).  This  is  the  North-East- 
Down  (NED)  convention. 

The  body  frame  (b- frame)  has  its  origin  at  the  center  of  mass  of  a  vehicle. 
The  x ,  y,  and  z  axes  are  aligned  with  the  roll  (0),  pitch  ( 6 ),  and  yaw  (0)  axes, 
respectively.  The  three  possible  translations  about  the  body  axes  (x,  y,  and  z)  and 
the  three  rotations  (roll,  pitch,  and  yaw)  are  the  six  degrees  of  freedom  for  the  aircraft. 
Figure  2.2  illustrates  a  typical  aircraft  b-frame.  Since  this  research  will  concern  two 
aircraft  flying  in  formation,  the  body  frames  of  each  will  be  designated  based  on  their 
formation  position.  Coordinates  relative  to  the  lead  aircraft  (tanker)  body  frame  will 
be  denoted  with  an  L  (e.g.,  pL).  Coordinates  relative  to  the  wing  aircraft  (receiver) 
body  frame  will  be  denoted  with  a  W  (e.g.,  p'v). 
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Figure  2.2:  Example  aircraft  body  frame  illustration 

2. 4  Coordinate  Transformations 

With  the  large  number  of  available  reference  frames,  a  vector  will  often  need 
to  be  resolved  from  one  frame  into  another  frame.  This  brings  rise  to  the  coordinate 
transformation,  which  is  often  accomplished  by  use  of  the  direction  cosine  matrix 
(DCM).  The  DCM  from  frame  a  to  frame  b  is  denoted  by  Cba.  As  shown  in  [17],  this 
3x3  matrix  is  the  product  of  every  unit  basis  vector  in  frame  a  with  every  unit  basis 
vector  in  frame  b.  Each  element  in  the  DCM,  CV,-,  is  the  cosine  of  the  angle  between 
the  i  axis  in  the  b  frame  with  the  j  axis  in  the  a  frame. 

A  vector  in  basis  b  can  easily  be  calculated  by  pre-multiplying  the  vector  of 
basis  a  by  the  DCM  as  shown  below: 


r  =  C 


bj.a 


(2.31) 
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2.5  WGS-84  Reference  System 

The  Earth  is  a  highly  irregular  body  and  various  methods  have  been  employed  to 
model  the  Earth’s  shape.  Though  similar  to  a  sphere,  the  Earth  bulges  at  the  equator 
and  is  more  flattened  at  the  poles.  For  precise  navigation,  modeling  the  Earth  as  a 
sphere  leads  to  errors.  The  ellipsoid  has  been  adopted  as  the  most  common  datum 
to  model  the  Earth’s  shape.  A  geodetic  datum  defines  a  three  dimensional  (3D) 
coordinate  system  to  identify  a  specific  location  in  terrestrial  space.  These  dimensions 
consist  of  geodetic  latitude  and  longitude  and  height  (above  the  ellipsoid)  [17]. 

In  [1],  The  World  Geodetic  System  Committee  created  a  datum  which  approx¬ 
imates  the  shape  of  the  Earth  over  the  entire  globe.  Created  in  1984,  this  datum  is 
known  as  WGS-84.  For  mapping,  charting,  geopositioning,  and  navigation,  WGS-84 
is  the  most  widely  used  global  geodetic  reference.  With  the  origin  at  the  center  of 
the  Earth,  the  z  axis  aligns  with  the  axis  of  rotation  and  the  x  and  y  axes  align  with 
the  equatorial  plane.  This  is  coincident  with  the  e-frame  described  in  the  previous 
section.  The  x  axis  points  toward  the  Zero  Meridian,  with  the  y  axis  orthogonal 
to  it.  From  Figure  2.1,  longitude  (l)  is  the  angle  between  the  local  meridian  plane 
and  the  Greenwich  meridian.  Latitude  (L)  is  the  angle,  in  the  local  meridian  plane, 
between  the  local  navigation  axes  and  the  Equator.  The  mathematical  model  which 
approximates  the  shape  of  the  Earth  is  known  as  the  ellipsoid.  Position  coordinates 
in  the  WGS-84  reference  frame  are  given  in  latitude,  longitude,  and  height  above  the 
ellipsoid  (HAE). 

2.6  Inertial  Navigation  Systems 

This  section  describes  the  basic  aspects  of  an  inertial  navigation  system  (INS). 
From  [21],  an  INS  is  composed  of  accelerometers  and  gyroscopes  to  sense  motion, 
as  well  as  a  computer.  With  an  initial  set  of  conditions  provided  by  the  user,  the 
INS  is  able  to  determine  motion  and  relative  position  based  on  integrations  of  the 
raw  measurements  from  the  devices.  These  raw  measurements  are  used  to  determine 
change  in  velocity,  change  in  attitude,  and  change  in  geographic  position. 
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2. 6. 1  Basic  Functions  of  Components.  INS  systems  are  standalone,  meaning 
they  do  not  require  input  from  any  outside  source.  Based  on  knowledge  of  initial 
conditions,  the  system  is  able  to  track  changes  in  velocity  and  position  without  the 
aid  of  an  external  reference.  There  are  two  different  physical  configurations  for  most 
INS  instruments,  gimbaled  and  strapdown.  Gimbaled  instruments  are  mounted  to  a 
platform  that  is  connected  to  gimbals,  which  allow  the  platform  to  stay  in  a  fixed 
orientation  independent  of  the  attitude  of  the  vehicle.  Based  on  inputs  from  the 
sensors,  torque  motors  rotate  the  platform  to  counteract  any  motion  by  the  vehicle. 
These  systems  are  highly  accurate,  though  they  require  a  high-level  of  maintenance 
due  to  complex  moving  parts.  Strapdown  systems,  just  as  the  name  suggests,  are  fixed 
to  the  body  frame  of  the  vehicle.  Typically,  one  gyroscope  and  one  accelerometer 
are  aligned  with  each  axis  to  measure  each  axis  of  rotation  and  each  translational 
acceleration.  Velocity  and  position  are  calculated  by  integrating  the  raw  data  from 
the  gyroscopes  and  accelerometers  [21], 

Gyroscopes  are  able  to  accurately  measure  the  rotational  accelerations  of  the 
vehicle  body  with  respect  to  inertial  space.  Accelerometers,  however,  measure  specific 
force,  which  is  the  sum  of  gravity  and  the  acceleration  of  the  vehicle  with  respect  to 
inertial  space.  Thus,  they  are  unable  to  separate  the  ever-present  effects  of  gravity. 
Hence,  an  INS  computer  must  subtract  the  local  gravitational  field  intensity  to  resolve 
true  acceleration  relative  to  inertial  space  [17]. 

The  effect  of  the  Earth’s  rotation  must  also  be  accounted  for  in  INS  position 
calculations.  The  Coriolis  equations  factor  in  the  effects  of  the  rotating  Earth  on  the 
moving  vehicle.  Though  not  discussed  in  detail  here,  they  are  an  important  factor  for 
INS  design  [17]. 

2.6.2  Mechanization.  Now  that  the  main  components  of  the  INS  section  are 
introduced,  the  mechanization  for  the  INS  can  be  discussed. 

As  discussed  in  [6],  the  INS  navigation  equations  are  used  to  describe  how 
the  position,  velocity,  and  attitude  are  updated  over  time.  Using  the  angular  rate 
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measurement,  ubb,  the  attitude  solution  is  updated  in  the  body  to  Earth-frame  DCM, 
Cb,  through  the  following  calculation 


C 


e 

b 


c  eA  -  nict 


(2.32) 


The  terms  f \\e  and  Qbb  are  the  skew  symmetric  matrices  of  the  Earth-rotation  vector 
and  the  IMU’s  angular-rate  measurement,  respectively. 

To  update  velocity,  one  must  note  that  the  centrifugal  and  Coriolis  components 
must  be  included  to  account  for  the  rotation  of  the  Earth-frame  axes  with  respect  to 
inertial  space.  The  velocity  update  is  computed  as: 


Kb  =  -  2  niv%  +  (2.33) 

The  acceleration,  a®b,  is  the  sum  of  both  the  measured  specific  force  (f4efe)  and  the 
gravitational  force  (7^).  Also,  one  can  define  the  term,  g|,  as  the  sum  of  the  centrifugal 
and  gravitational  accelerations: 


g b  ~  lib  ^ie^iePeb 


Thus,  making  these  two  substitutions  to  Equation  2.33  yields: 


Kb  =  f ib  +  gb(Peft)  -  2 neiev[ 


The  position  change  is  simply  equal  to  the  velocity: 


(2.34) 


(2.35) 


Pefo  =  Vefe 


(2.36) 
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2.1  Optical  Sensors 

This  section  describes  some  fundamental  concepts  regarding  camera  models.  For 
this  research,  a  pinhole  camera  model  was  assumed  for  the  digital  camera  onboard 
the  wing  (receiver)  aircraft.  This  is  consistent  with  the  flight  data  used  for  image 
processing  analysis,  which  will  be  discussed  in  detail  in  Chapters  III  and  IV. 

Figure  2.3  shows  a  typical  pinhole  camera  model.  The  image  plane  is  projected 
one  focal  length  (/)  in  front  of  the  lens.  As  shown  in  the  figure,  the  coordinate 
frame  for  the  camera  is  right-handed,  with  the  x  and  y  axes  pointing  down  and 
right,  respectively.  The  0  axis  points  perpendicular  to  the  camera  lens.  For  this 
research,  there  is  no  need  to  track  specific  pixels  in  any  images  (i.e. ,  the  location  of 
the  projection  on  the  image  plane). 


Sc 


Scene 


Figure  2.3:  Simple  Camera  Projection  Model.  The  virtual 

image  plane  is  placed  in  between  the  focal  point  and  the  scene, 
thus  the  image  is  not  inverted  as  typical  in  optical  models. 
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2.8  Image  Processing  Techniques 

During  a  refueling  operation,  the  digital  camera  onboard  the  UAS  will  take 
periodic  images  of  the  tanker  aircraft.  The  relative  navigation  system  will  generate 
predicted  images  of  the  tanker  based  on  the  current  state  estimate  and  covariance. 
These  predicted  images  will  be  compared  to  the  measured  image,  to  provide  a  virtual 
measurement  as  feedback  to  the  system.  More  details  of  how  the  system  creates  a 
predicted  image  will  be  discussed  in  Section  2.9 

For  analyzing  image  processing  techniques,  two  specific  formation  positions  were 
examined:  the  pre-contact  and  contact  positions.  At  the  pre-contact  position,  the  re¬ 
ceiver  aircraft  is  approximately  50  feet  behind  and  10  feet  below  the  tip  of  the  refueling 
boom.  In  the  contact  position,  the  receiver  aircraft  is  positioned  just  below  the  fixed 
boom,  with  its  refueling  receptacle  directly  in  line  with  the  boom  extension.  This 
is  the  optimal  refueling  location.  Figure  2.4  presents  all  typical  formation  positions 
during  aerial  refueling. 


Tanker  Body  Y  Position,  ft 


Figure  2.4:  Diagram  of  typical  formation  positions  during 

aerial  refueling.  This  research  will  focus  on  the  pre-contact  and 
contact  positions  [3]. 
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Figure  2.5  shows  example  infrared  images  from  the  receiver  camera  perspective 
while  in  the  pre-contact  and  contact  positions.  The  relative  positioning  accuracy 
while  the  receiver  aircraft  is  in  pre-contact,  contact,  and  transitioning  between  these 
two  points  is  absolutely  critical.  Due  to  the  close  proximity  of  the  two  aircraft  dur¬ 
ing  this  portion  of  the  refueling  operation,  any  positioning  errors  could  be  a  safety 
hazard.  Thus,  these  are  optimal  positions  to  evaluate  the  accuracy  of  an  image-aided 
algorithm.  The  simulations  presented  in  Chapter  IV  will  focus  on  these  formation 
positions. 

The  key  goal  for  image  processing  in  this  research  is  to  determine  some  type  of 
metric  between  the  predicted  and  true  images.  The  system  must  have  some  way  to 
determine  if  there  are  errors  between  the  two  images,  and  how  large  those  errors  are. 
Exactly  how  the  system  uses  these  errors  as  measurement  updates  will  be  discussed 
in  Section  3.2.2.  Three  image  comparison  metrics  are  considered  in  this  research: 
a  sum-squared  difference  of  intensity,  a  magnitude  of  gradient  comparison,  and  a 
gradient  with  threshold  technique.  These  approaches  attempt  to  quantify  the  error 
between  a  given  predicted  image  and  the  truth  image.  Each  approach  is  described  in 
the  paragraphs  that  follow. 


(a)  (b) 

Figure  2.5:  (a)  Example  infrared  image  of  the  tanker  as  seen 

from  the  receiver  at  the  pre-contact  position. 

(b)  Example  infrared  image  of  the  tanker  as  seen  from  the  re¬ 
ceiver  at  the  contact  position 
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2.8.1  Sum  Squared  Difference.  The  two  images  are  first  converted  to 
grayscale  and  the  respective  image  intensities  are  normalized.  This  normalization 
is  computed  by  first  subtracting  the  mean  image  value  (I)  from  each  pixel,  then  di¬ 
viding  by  the  standard  deviation  (a)  of  the  pixel  intensities.  For  a  given  image,  I.  the 
normalization  calculation  is  shown  in  Equation  2.37. 


I 


norm 


i-l 

(7 1 
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In  order  to  avoid  aliasing  effects,  the  image  is  passed  through  a  filter  that  applies  a 
Gaussian  blur  to  smooth  the  edges  of  the  aircraft.  This  Gaussian  blur  (mean  =  0, 
a  =  1  pixel)  will  be  applied  during  all  image-processing  techniques. 

Now,  to  find  the  sum-squared  difference,  c,  between  an  image  I  and  another 
image  N  (each  with  n  pixels),  the  following  expression  is  calculated: 

n 

c  =  J](I i  -  N,;)2  (2.38) 

i=  1 


Conceptually,  the  sum-squared  difference  approach  simply  compares  the  differ¬ 
ences  in  intensity  between  the  pixels  in  one  image  with  the  corresponding  pixels  in 
another  image.  This  technique  would  work  well  when  there  is  a  distinct  contrast  in 
intensity  between  the  aircraft  portion  of  the  image  and  the  background  portion  of  the 
image.  This  is  indeed  true  in  the  case  of  infrared  imaging  and  the  images  generated 
by  predictive  rendering  (Section  2.9). 


2.8.2  Magnitude  of  Gradient.  A  second  approach  considered  for  image 
correlation  is  to  convert  the  images  to  gradient  space.  A  gradient  is  a  mathemati¬ 
cal  calculation  to  determine  the  change  in  image  intensity  between  adjacent  pixels. 
Areas  containing  a  sudden  change  in  pixel  intensity  indicate  edges  or  corners.  In 
this  research,  these  typically  exist  at  the  edges  of  the  tanker  aircraft  in  the  image. 
Applying  a  special  image  processing  filter  in  Matlab®  known  as  the  Prewitt  filter 
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(Equation  2.39)  to  the  image,  one  can  determine  the  horizontal  gradient  matrix  (Gx) 
for  the  image. 
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(2.39) 


Simply  taking  the  transpose  of  this  filter  will  yield  the  gradient  matrix  for  the  vertical 
direction  (Gy).  The  overall  gradient  magnitude  matrix  in  the  image  is  calculated  as: 
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(2.40) 


This  magnitude  is  then  normalized  across  the  image,  similar  to  the  calculation  above 
for  sum-squared  difference. 
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2.8.3  Magnitude  of  Gradient  with  Threshold.  As  a  third  approach,  the 
same  magnitude  of  gradient  is  computed  for  each  image,  but  the  magnitudes  are 
then  filtered  based  on  a  threshold  technique.  Through  experimentation,  a  threshold 
is  found  that  will  amplify  strong  areas  of  gradient  magnitude.  Each  pixel  in  the 
image  is  compared  to  the  threshold  value  and  is  set  to  a  value  of  zero  if  below  the 
threshold,  or  set  to  a  value  of  one  if  above  the  threshold.  This  will  remove  the  few 
small,  unwanted  edges  (any  background  clutter  or  small  features  on  the  surface  of  the 
aircraft),  as  well  as  amplify  the  desired  edges  of  the  aircraft  outline.  A  sample  image 
(after  implementing  a  threshold)  is  shown  in  Figure  2.6. 
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Figure  2.6:  Example  image  after  applying  magnitude  of  gra¬ 
dient  with  threshold.  The  threshold  on  the  gradient  amplifies 
the  edges  of  the  tanker  while  removing  any  unwanted  gradient 
due  to  background  clutter  or  lens  reflection. 

2.9  Predictive  Rendering  Algorithm 

2. 9. 1  Motivation.  The  basic  premise  of  vision-aided  navigation  is  to  extract 
as  much  information  as  possible  from  digital  images.  In  this  research  application, 
the  relative  distance  and  attitude  between  the  two  aircraft  must  be  determined  based 
on  images  taken  from  the  wing  aircraft  imaging  sensor.  In  order  to  extract  this 
information,  there  must  be  some  initial  prediction  of  what  the  image  should  look  like. 
As  discussed  in  the  previous  section,  an  image  processing  algorithm  will  compare 
a  predicted  image  with  the  truth  image  in  order  to  provide  feedback  to  the  filter. 
Without  a  prediction  for  comparison,  the  truth  image  alone  does  not  reveal  any 
specific  metric  of  relative  distance  and  attitude. 

The  predictive  rendering  algorithm  generates  a  predicted  image  using  a  three- 
dimensional  structural  model  of  the  lead  aircraft  combined  with  an  a  priori  estimate 
of  the  relative  pose  between  the  imaging  sensor  and  the  lead  aircraft. 
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Many  image-aided  navigation  techniques  for  aerial  refueling  focus  on  predicting 
and  tracking  specific  features  (edges  or  corners)  on  the  tanker  aircraft.  This  can  be  a 
computational  burden  on  the  system  because  of  the  required  processing  to  track  and 
associate  multiple  features  throughout  the  filter  execution.  Furthermore,  there  is  the 
possibility  of  misassociating  a  detected  feature  with  an  incorrect  predicted  feature. 
If  the  pose  or  range  estimate  is  significantly  wrong,  the  portion  (or  window)  of  the 
image  where  the  filter  predicts  a  specific  feature  could  potentially  contain  a  different 
feature  of  the  aircraft.  This  can  lead  to  even  larger  errors  in  the  relative  navigation 
solution  if  the  feature  detection  algorithm  associates  the  wrong  detected  feature  with 
a  predicted  feature.  In  other  words,  potential  errors  exist  with  predicting  a  certain 
feature  in  a  small  window  of  the  image  and  finding  the  wrong  feature.  Similarly, 
predicting  a  feature  in  a  specific  portion  of  the  image  and  detecting  no  features  also 
leads  to  errors  [15]. 

In  order  to  simplify  this  problem,  the  predictive  rendering  approach  seeks  to 
generate  the  entire  image  of  the  tanker  and  compare  it  with  the  entire  image  taken 
from  a  digital  camera  using  image-processing  techniques.  This  eliminates  the  need 
to  track  specific  features  on  the  aircraft  and  also  eliminates  the  possibility  of  feature 
misassociations.  Additionally,  this  technique  does  not  require  modification  of  the 
tracking  processes  when  certain  features  of  the  tanker  are  no  longer  visible  to  the 
camera.  For  example,  when  transitioning  from  the  pre-contact  to  the  contact  position, 
the  outer  extremities  of  the  lead  aircraft  (i.e.,  wingtips  and  tail  section)  slowly  move 
out  of  the  camera’s  field  of  view.  A  feature  detection  algorithm  would  have  to  discard 
these  features  as  they  moved  out  of  the  image  and  rely  on  new  features  (if  they  exist 
and  are  detectable).  Images  of  the  tanker  taken  during  these  two  formation  positions 
are  shown  in  Figure  2.5,  presented  earlier  in  this  chapter. 

Additionally,  the  predictive  rendering  technique  will  give  the  relative  navigation 
algorithm  a  modular  design.  The  fundamental  algorithm  will  stay  the  same  regardless 
of  which  lead  aircraft  the  wing  aircraft  is  following  in  the  formation.  A  3D  model 
of  a  different  tanker  aircraft  (e.g.  KC-10  Extender)  could  be  substituted  in  for  the 
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KC-135R,  and  the  relative  navigation  solution  would  function  the  same.  Furthermore, 
as  long  as  a  3D  model  is  available,  virtually  any  aircraft  could  be  used  as  the  lead 
aircraft  for  other  applications  of  autonomous  formation  flight. 

2.9.2  Implementation.  In  this  research,  the  rendering  algorithm  was  imple¬ 
mented  in  Matlab®  using  the  Virtual  Reality  (VR)  Toolbox  and  a  three  dimensional 
model  of  the  KC-135R  Stratotanker.  The  VR  toolbox  uses  OpenGL  to  render  scenes 
of  interest.  The  VR  environment  simply  needs  a  translation  vector  and  rotation  vec¬ 
tor  to  adjust  the  location  of  the  3D  model  within  the  window.  This  approach  was 
used  clue  to  the  ease  with  which  Matlab®  simulations  could  be  developed  that  can 
directly  interface  with  the  rendering  tool.  Though  the  VR  toolbox  cannot  render  at 
a  frame  rate  as  high  as  more  powerful,  commercial  rendering  tools,  it  has  adequate 
performance  to  prove  the  concept  of  this  research.  It  is  not  an  adequate  approach  for 
real-time  operation.  For  example,  if  the  system  captured  a  truth  image  at  a  rate  of 
1  Hz,  the  rendering  tool  would  be  required  to  generate  many  predicted  images  at  a 
very  high  frame  rate.  These  images  would  then  have  to  be  processed  and  used  for  an 
update  to  the  filter,  all  in  enough  time  for  the  system  to  receive  another  image  update 
one  second  later.  More  on  this  update  process  will  be  discussed  in  Section  3.2.2. 

The  relative  pose  of  the  lead  aircraft  is  defined  as  the  position  and  orientation  of 
the  lead  body  frame  relative  to  the  camera  frame  on  the  wing  aircraft.  This  position 
and  orientation  is  calculated  by: 

s'  =  c;,,  [c"c;(p|  -  Pew)  -  pf]  (2.42) 

Cl  =  C'b,C«'C2  (2.43) 

These  relative  position  and  orientation  parameters  are  then  provided  to  the 
predictive  rendering  algorithm.  The  VR  tool  then  generates  a  3D  rendering  of  the 
predicted  image  for  any  desired  translation  and  rotation.  A  sample  rendering  is  shown 
in  Figure  2.7. 
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Figure  2.7:  Sample  3D  Rendering  of  Tanker  Aircraft. 

Mat  lab®  Virtual  Reality  environment  renders  a  predicted  image 
of  the  lead  tanker  from  the  perspective  of  the  receiver  aircraft. 

The  simulation  environment  was  validated  using  data  collected  during  flight 
tests  by  the  Northrop  Grumman  Corporation.  These  data  sets  included  time-space 
positioning  information  (TSPI)  data  for  both  the  receiver  and  tanker  aircraft,  as  well 
as  infrared  imaging  from  a  camera  mounted  inside  the  nosecone  of  the  receiver.  The 
camera  was  mounted  with  a  pitch  angle  of  +28°  in  order  to  keep  the  tanker  in  its 
held  of  view  while  the  receiver  aircraft  was  below  the  tanker  in  the  contact  position. 
The  specific  type  of  imaging  sensor  used  was  Long  Wave  Infrared  (LWIR).  This  type 
was  chosen  during  these  data  collection  flights  due  to  the  ease  of  locating  the  tanker 
aircraft  in  the  infrared  image  in  all  lighting  and  weather  conditions  [3].  An  image  of 
the  camera  mounted  in  the  nose  of  the  receiver  aircraft  is  shown  in  Figure  2.8. 

A  sample  predicted  image  from  the  simulation  is  shown  in  Figure  2.9.  Here, 
the  predicted  location  of  the  tanker  (with  an  intentional  error  for  easier  distinction) 
is  rendered  in  the  same  figure  as  the  true  infrared  image  of  the  tanker.  The  predicted 
and  measured  images  must  be  compared  in  order  to  determine  the  accuracy  of  the 
a  priori  state  estimate.  This  comparison  metric  is  generally  termed  as  the  residual 
between  predicted  and  measured  values  and  is  described  in  Section  3.2.2. 
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Figure  2.8:  Orientation  of  the  LWIR  Camera  when  mounted 
in  the  nose  of  the  receiver  aircraft  during  data  collection  flights. 
The  camera  has  a  pitch  angle  of  +28°. 
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Figure  2.9:  Rendered  Image  Prediction  (Semi-Transparent) 

with  True  Infrared  Image.  A  semi-transparent  rendering  of  the 
predicted  tanker  location  is  superimposed  on  an  actual  infrared 
image  from  the  receiver  camera. 


26 


2.10  Previous  Research  Efforts 

This  section  highlights  previous  work  done  regarding  autonomous  aerial  refuel¬ 
ing,  including  all  significant  approaches  (GPS,  vision-aided,  etc.)  to  date. 

2.10.1  Ross/Spinelli  Formation  Flight  Controller.  The  goal  of  the  Ross  and 
Spinelli  research  was  to  develop  a  system  that  would  demonstrate  fully-autonomous 
formation  flight  and  accurate  relative  position  to  a  lead  tanker  aircraft.  This  involved 
maintaining  proper  position  through  all  major  formations  associated  with  inflight 
refueling  (wing,  contact,  and  pre-contact)  [13, 16].  The  fundamental  approach  for 
positioning  (developed  by  Spinelli)  was  to  use  differential  GPS  (DGPS)  between  the 
two  aircraft.  Using  a  wireless  datalink,  the  receiver  used  the  tanker  position  as  a  local 
truth  reference,  and  through  double-differencing,  obtained  its  position  relative  to  the 
tanker  at  centimeter  level  accuracy  [16]. 

Ross  successfully  created  a  flight  controller  that  maintained  relative  position 
optimized  for  15°  of  bank,  but  capable  of  up  to  30°.  The  system  under  test  was  a 
Calspan  Learjet  as  the  receiver  aircraft  with  a  C-12  aircraft  as  the  surrogate  tanker. 
Figure  2.10  shows  the  two  aircraft  in  formation  during  one  of  the  flight  tests.  This 
research  successfully  demonstrated  that  GPS  is  one  valid  approach  for  accurate  rel¬ 
ative  positioning.  Ross  also  noted  that  if  an  optical  system  could  keep  the  receiver 
aircraft  within  a  five  foot  envelope  of  the  nominal  refueling  position,  the  human  boom 
operator  can  correct  the  error  and  still  successfully  make  contact  with  the  refueling 
boom  [13]. 

2.10.2  Spencer  Optical  Tracking  for  Relative  Positioning.  Using  feature 
extraction,  matching,  and  validation,  Spencer  aimed  to  determine  relative  position 
between  two  aircraft  using  optical  techniques.  The  fundamental  feature  extraction 
technique  was  the  Harris  corner  detector.  For  navigation,  the  system  attempts  to 
find  and  propagate  the  most  likely  association  feature  to  the  next  frame  using  a 
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Figure  2.10:  Image  of  Calspan  Learjet  and  C-12  surrogate 

tanker  aircraft  during  Ross/Spinelli  Formation  Flight  Test 

traditional  Kalman  Filter.  The  system  then  tries  to  match  a  detected  corner  or  edge 
from  the  image  with  a  predicted  feature  of  the  target  [15]. 

The  system  used  pose  estimates  from  the  inertial  measurement  unit  (IMU) 
coupled  with  the  measurement  update  from  the  Kalman  filter’s  predictions  of  the 
feature  locations.  The  radial  projection  error  (between  projected  features  and  actual 
locations)  was  approximately  five  pixels  for  relative  distance  between  60-300  feet, 
and  approximately  8-12  pixels  inside  of  60  feet  during  Spencer’s  flight  tests.  Due 
to  this  error,  a  different  type  of  feature  extraction  or  image  matching  method  could 
potentially  allow  this  approach  to  more  accurately  determine  the  relative  position  and 
also  be  more  computationally  tractable  [15]. 

2.10.3  Vision-Based  Sensor  and  Navigation  System.  Another  visual  track¬ 
ing  system  for  refueling  was  developed,  in  concept,  at  Texas  A&M  University,  called 
the  Vision-Based  Sensor  and  Navigation  System,  or  VizNav.  The  approach  here  is  to 
use  the  probe  and  drogue  technique  of  inflight  refueling,  with  several  infra-red  light 
emitting  diodes  (LEDs)  affixed  to  the  drogue  chute  in  a  known  orientation.  When 
the  receiver  aircraft  approaches  the  drogue,  the  system  aboard  the  receiver  aircraft 


establishes  a  wireless  datalink  with  the  electronics  mounted  on  the  drogue.  The  sys¬ 
tem  then  sends  commands  to  the  drogue  electronics  which  then  activate  each  LED 
in  a  precise,  pre-determined  sequence.  The  optical  sensor  aboard  the  receiver  air¬ 
craft  can  detect  and  distinguish  each  LED  and  thus  determine  a  relative  distance  and 
orientation  to  the  drogue  [18]. 

Though  a  very  interesting  and  seemingly  robust  approach,  this  system  does  de¬ 
mand  the  ability  to  establish  a  reliable  datalink  between  the  aircraft  and  the  drogue 
system.  A  suitable  corollary  with  boom  refueling  isn’t  readily  apparent,  since  with 
the  boom  approach  the  receiver  aircraft  does  not  have  a  visual  target  (the  drogue)  to 
intercept.  Furthermore,  this  approach  would  require  undesirable  and  costly  modifica¬ 
tions  (installing  LEDs)  to  the  tanker  fleet  [18]. 

2.10.4  DARPA/NASA  Dry  den  Probe  and  Drogue  Refueling.  The  autonomous 
aerial  refueling  system  developed  under  a  joint  Defense  Advanced  Research  Projects 
Agency  (DARPA)  and  National  Aeronautics  and  Space  Administration  (NASA)  test 
at  the  NASA  Dryden  test  facility  involved  GPS,  vision,  and  inertial  sensors  to  de¬ 
termine  relative  position  between  a  refueling  drogue  from  a  Boeing  707/300  and  an 
F/A-18  aircraft. 

The  system  primarily  calculated  the  position  of  the  system  using  the  combined 
relative  GPS  and  INS  measurements.  Once  within  visual  reference  (approximately  80- 
120  ft),  the  vision  system  would  track  the  drogue  using  an  a  priori  model  of  the  drogue 
shape.  The  test  established  an  autonomous,  hands  off,  repeatable  demonstration  of 
positive  contact  between  the  receiver  aircraft  and  refueling  receptacle  [7]. 

2.10.5  Perugia  University  and  West  Virginia  University.  Perugia  University 
and  West  Virginia  LIniversity  have  many  research  contributions  in  the  autonomous 
aerial  refueling  problem  [4,5,14],  These  efforts  focus  on  comparing  different  machine 
vision  algorithms  for  feature  detection  on  the  lead  tanker  aircraft  and  how  to  develop 
robust  controllers  for  autonomously  controlling  the  receiver  aircraft  in  the  formation. 
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These  vision-aided  approaches  all  consider  tracking  specific  features  on  the  lead  air¬ 
craft,  sometimes  even  using  infrared  optical  markers  in  a  known  orientation  on  the 
tanker.  The  simulations  from  Perugia  and  West  Virginia  Universities  show  favorable 
results,  but  the  research  presented  in  this  thesis  seeks  to  avoid  tracking  specific  fea¬ 
tures  of  the  tanker  and/or  rely  on  modifying  the  tanker  fleet  by  installing  optical 
markers. 

2.11  Summary 

This  chapter  has  highlighted  key  mathematical  and  navigation  background  that 
will  be  necessary  for  thorough  understanding  of  Chapter  111.  The  EKF  basic  discussion 
will  be  used  to  develop  relative  navigation  filters  that  use  INS  sensors  for  propagation 
and  images  as  virtual  measurement  updates. 

The  related  research  presented  in  this  chapter  summarizes  key  work  and  suc¬ 
cesses  already  accomplished  to  solve  the  autonomous  aerial  refueling  problem.  The 
algorithm  presented  in  the  next  chapter  seeks  to  complement  these  efforts  to  obtain 
a  robust,  image- aided  navigation  solution. 


30 


III.  Algorithm  Development  and  Implementation 

This  chapter  will  describe  the  algorithm  used  to  develop  a  vision-aided  relative 
navigation  system.  First,  the  dynamics  model  is  developed  for  both  aircraft  in 
the  formation.  Several  image  processing  techniques  are  presented  to  determine  the 
accuracy  of  the  rendering.  Finally,  the  Extended  Kalman  Filter  implementation  is 
discussed,  including  the  algorithm  to  use  images  as  the  measurement  update. 

3.1  Modeling  Aircraft  Dynamics 

Before  constructing  a  filter  to  estimate  the  relative  navigation  solution,  the 
system  dynamics  must  be  accurately  modeled,  ft  is  assumed  that  an  INS  is  available 
on  the  wing  aircraft,  so  the  wing  dynamics  are  modeled  based  on  the  nonlinear  INS 
mechanization  equations  and  associated  statistical  error  model. 

3.1.1  Wing  Aircraft  Dynamics.  The  model  used  in  this  research  builds  upon 
the  mechanization  development  presented  in  [19].  The  accelerometer  measurement 
model  is: 

{w  =  fw  +  aw  +  ww  (3.!) 

Recall  from  section  2.3  that  the  L  and  W  notation  refers  to  the  lead  and  wing  body 
frames,  respectively.  The  measured  specific  forces  in  the  wing  frame  (f,1^ )  are  the  sum 
of  the  true  specific  forces  in  the  wing  frame  ( fw ),  the  accelerometer  biases  (aiv),  and 
additive  white  Gaussian  noise  (w^  ).  The  gyroscope  measurement  model  is: 

=  a#  +  b"'  +  wf  (3.2) 

The  measured  body  rotation  rates  with  respect  to  inertial  space  are  the  sum 

of  the  true  body  rotation  rates  (u>^),  the  gyroscope  biases  (bw ),  and  additive  white 
Gaussian  noise  (w xbv).  The  biases  for  both  the  accelerometers  and  the  gyroscopes  are 
modeled  as  first-order  Gauss-Markov  processes: 

a-  =  __V'+w»;M  (3.3) 

•  n. 
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(3.4) 


l'W  1  u  W  I  w 

b  = - b  +  Wu 

n 

The  time  constants  for  the  accelerometer  and  gyroscope  are  represented  by  ra  and  Tb, 
respectively. 

Fully  developed  in  [19],  the  error-state  model  is  presented  in  the  following  equa¬ 
tions.  First,  the  attitude  error  states  are  defined  in  terms  of  small  angles  about  the 
north,  east,  and  down  directions  of  the  local  navigation  frame  as: 

fpn 

(3-5) 

The  final  attitude  error  equations  (in  vector  form  from  [19])  are  determined  to  be: 

ip  =  —  (C  gUjfe)x  ip  —  C^/b11  —  (3-6) 

The  x  symbol  following  a  vector  indicates  the  skew  symmetric  form.  The  velocity 
and  position  error  equations  (in  vector  form  from  [19])  are  determined  to  be: 

8vn  =  CneGCeJpn  -  2CneneieCeJvn  +  (P x)?/>  +  C”a6  +  Cnbwba  (3.7) 

hp«  =  8vn  (3.8) 

3.1.2  Lead  Aircraft  Dynamics.  If  assuming  that  the  lead  INS  measurements 
are  available  for  the  navigation  solution,  the  dynamics  are  modeled  using  the  same 
error  equations  as  presented  above  in  the  wing  dynamics  section.  Since  this  research 
will  focus  on  a  completely  passive  system,  the  INS  measurements  are  not  available 
and  the  dynamics  are  modeled  using  a  stochastic,  kinematic  motion  model.  With 
the  knowledge  that  the  lead  flight  behavior  will  be  relatively  benign  during  a  close- 
proximity  refueling  maneuver,  several  of  the  states  can  be  modeled  as  first-order 
Gauss-Markov  processes  driven  by  WGN. 


32 


The  x  velocity  dynamics  in  the  lead  body  frame  are  modeled  as: 


+  <  (3-9) 

Tv 

The  variables  are  tracked  in  the  lead  body  frame  in  order  to  make  the  first-order 
Gauss-Markov  process  a  valid  modeling  tool.  The  y  and  z  velocities  in  the  lead  body 
frame  should  have  very  minute  variations  and  each  can  simply  be  modeled  by  a  WGN 
source. 

i’y  =  '.'t,  (3.10) 

y  =  <  (3.11) 

Similar  to  the  x  velocity,  the  roll  ((f>L)  and  pitch  (0L)  state  dynamics  are  also  modeled 
as  first-order  Gauss-Markov  processes. 


<Pl  = - (Pl  +  wh 

t4>l 


(3.12) 


1 

t0l 


0l  = - 0L  +  WeT 


(3.13) 


The  time  constants  for  the  roll  and  pitch  rates  are  and  TgL,  respectively. 

The  heading  (ipL,)  differential  equation  is  a  function  of  the  amount  of  roll,  the 
gravity  vector  (g),  and  the  x  velocity  as  follows: 


= 


g  tan  cpL 


(3.14) 


JXL 


Now  to  convert  to  the  error  state  form,  mixed  partial  derivatives  of  each  equa¬ 
tion  above  are  taken  with  respect  to  each  of  the  state  variables  (commonly  called 
the  Jacobian).  The  lead  velocity  error  state  variable  dynamics  (in  vector  form)  are 
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computed  as: 


— 0  0 

Tvx 

0  0  0 

0  0  0 


(3.15) 


The  lead  position  error  state  variable  dynamics  (in  the  navigation  frame)  are  simply 
a  function  of  the  velocity: 


6pl  =  CnL8vLL 

The  attitude  error  state  variable  dynamics  are  as  follows: 
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(3.17) 


3.1.3  State-Space  Model.  With  all  of  the  error  state  differential  equations 
now  written  for  the  entire  system,  the  overall  state-space  model  can  be  formed.  This 
will  be  the  nonlinear  model  used  in  the  propagation  step  of  the  Kalman  filter. 

The  error  state  vector,  <5x,  consists  of  the  position  (<5p(^),  velocity  (5v(^),  atti¬ 
tude  accelerometer  bias  (<5aM/),  and  gyroscope  bias  (hbu  )  for  the  wing  vehicle, 

as  well  as  the  position  (dpi),  velocity  (<5v£),  and  attitude  (i/?L  =  [84>l,  80l,  8'iPl]T) 
for  the  lead  vehicle.  This  yields  a  twenty-four  element  vector  of  the  form: 


5x  = 


dPw 

■*/v 

Sbw 

<*P  L 

K 


(3.18) 
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The  noise  term  consists  of  the  WGN  sources  that  drive  the  accelerometer  measurement 


(w^),  the  gyroscope  measurement  (w^),  the  accelerometer  bias  (w«Was)>  and  the 
gyroscope  bias  (w^a  )  all  for  the  wing  aircraft,  as  well  as  the  velocity  noise  (w %el), 
and  attitude  noise  (w^t  =  [w^wg,  0]T)  for  the  lead  aircraft.  These  elements  are 


combined  to  form: 


(3.19) 
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The  overall,  augmented  state-space  error  model  can  then  be  written  as  follows: 
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using  the  following  substitutions: 
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3.2  Image- Aided  Filter  Implementation 

This  section  will  highlight  the  key  details  of  how  the  filter  propagates  and  up¬ 
dates  during  operation. 

3.2.1  Extended  Kalman  Filter  Propagation.  As  discussed  in  Section  2.2.2, 
the  Extended  Kalman  Filter  propagates  from  one  measurement  update  to  the  next 
by  linearizing  the  trajectory  about  the  current,  nominal  state  estimate.  After  a  mea¬ 
surement  update  is  taken,  the  new  state  estimate  is  used  to  linearize  the  dynamics 
presented  in  the  previous  section.  This  linear  relationship  is  then  propagated  to  the 
next  time  epoch,  to  yield  the  a  priori  estimate.  The  filter  is  then  ready  to  receive  an 
update  based  on  the  image  captured  from  the  digital  camera. 

3.2.2  Image  Update.  Section  2.8  considered  three  methods  that  determine 
an  error  between  the  true  image  collected  by  the  digital  camera,  and  the  predicted 
image  generated  from  the  a  priori  state  estimate.  This  error  value  is  simply  a  nu¬ 
meric  number,  ft  quantifies  how  well  the  two  images  match,  but  has  no  information 
regarding  which  degrees  of  freedom  have  significant  errors.  Using  only  the  error,  one 
can  not  readily  deduce  which  direction(s)  to  move  the  optimal  state  estimate  so  that 
it  will  yield  a  prediction  that  more  closely  matches  the  true  image. 

By  perturbing  the  estimate  (in  all  six  degrees  of  freedom)  about  the  nominal 
and  generating  multiple  images,  the  image  update  can  determine  in  which  direction(s) 
the  estimate  differs  from  the  truth.  Before  these  perturbations  are  started,  there  must 
be  some  metric  to  regulate  how  large  the  perturbations  in  each  degree  of  freedom  will 
be.  If  the  initial  estimate  is  significantly  wrong,  then  small  perturbations  might  not 
give  the  algorithm  observability  into  the  appropriate  directions  that  will  correct  the 
estimate.  Also,  if  the  a  priori  estimate  is  very  accurate,  then  perturbations  that  are 
very  large  will  not  provide  observability  into  any  very  small  error  corrections  that  are 
needed. 
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The  state  error  covariance  matrix  (Pxx)  is  used  to  determine  how  large  the 
perturbation  values  should  be.  States  that  have  large  uncertainty  should  obviously 
have  larger  perturbations  than  those  with  smaller  uncertainty.  By  taking  the  Cholesky 
decomposition  of  this  matrix,  the  perturbation  values  for  each  state  can  be  found  that 
include  any  cross-coupling  covariance  effects  between  states.  Mat  lab®  has  a  built-in 
function  (. cholQ )  to  easily  yield  the  Cholesky  decomposition: 

chol(Pxx)T  =  y/¥7x  (3.22) 


Each  degree  of  freedom  is  then  perturbed  about  the  nominal  state  estimate,  using  the 
magnitudes  found  from  the  Cholesky  decomposition.  These  magnitudes  are  applied 
as  both  positive  and  negative  values,  so  that  the  estimate  is  perturbed  in  all  possible 
directions.  Mathematically,  each  perturbation  state  vector  (x*)  is  a  linear  function  of 
each  row  of  the  Cholesky  decomposition  plus  the  nominal,  a  priori  state  estimate  (x): 


x 


i 


+  X 


(3.23) 


After  each  individual  perturbation  generates  a  new  predictive  rendering,  an  image 
is  captured  and  compared  with  the  nominal,  a  priori  image.  This  yields  an  error 
gradient  in  all  directions  about  the  nominal,  serving  as  the  observability  function  (or 
measurement)  of  each  perturbation: 


z  i  =  h[±i]  (3.24) 

Once  every  gradient  is  found  (for  both  positive  and  negative  perturbations),  the  values 
are  compared  to  find  the  direction  that  produces  the  gradient  value  that  is  the  lowest 
(i.e.,  most  negative).  In  other  words,  moving  the  predicted  image  in  this  direction 
results  in  the  most  error  reduction  of  any  of  the  perturbations.  The  state  estimate 
is  adjusted  by  the  amount  of  the  best  perturbation,  and  the  process  repeats.  The 
estimate  is  then  perturbed  about  this  new  nominal  value  to  again  find  the  direction(s) 
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that  yields  the  greatest  reduction  in  error.  The  algorithm  tracks  the  translation  and 
rotation  adjustments  from  the  a  priori  estimate  throughout  the  process. 

As  this  process  cycles,  the  predicted  image  should  gradually  and  methodically 
move  so  that  it  closely  aligns  with  the  true  image  from  the  digital  camera.  Once 
the  predicted  and  true  images  are  closely  aligned,  then  the  normal  perturbations 
will  not  yield  a  direction  in  which  the  error  will  decrease.  Moving  the  prediction  in 
any  direction  will  only  increase  the  error.  At  this  point,  the  perturbation  values  are 
reduced  to  a  fraction  (10%)  of  their  original  value,  and  the  process  above  is  repeated. 
This  is  a  method  to  fine  tune  the  estimate  and  lock  in  the  predicted  image  with  the 
true  image  as  closely  as  possible. 

After  the  fine  tuned  perturbations  reach  the  point  where  any  movement  yields 
an  increase  in  error,  then  the  prediction  should  be  closely  correlated  with  the  true 
image.  This  location  is  indeed  a  minimum,  but  it  is  potentially  only  a  local  minimum. 
It  is  possible  that  the  global  minimum  could  be  at  a  different  position.  Using  the  state 
error  covariance  as  the  scale  with  which  to  perturb  helps  to  mitigate  this  possibility. 
If  the  covariance  values  are  high,  then  the  perturbations  should  be  high  enough  to  find 
the  global  minimum.  If  the  covariance  is  low,  indicating  the  nominal  estimate  is  near 
the  truth,  then  small  perturbations  should  find  the  global  minimum.  The  cumulative 
translation  and  rotation  adjustments  from  the  a  priori  estimate  to  the  closest  match 
yield  the  measured  location  and  orientation  of  the  lead  aircraft. 

As  shown  in  Equations  2.13  and  2.18  (for  both  the  traditional  and  Extended 
Kalman  Filters)  a  measurement  uncertainty  matrix,  R,  must  be  included  as  part  of 
the  measurement  update  step.  In  order  to  find  the  uncertainty  of  the  measurements 
taken  from  the  image  perturbations,  the  shape  of  the  error  curves  must  be  known. 
Every  error  value  is  stored  as  the  algorithm  perturbs  the  estimate  about  the  nominal 
until  finally  locking  on  to  the  best  matching  image.  These  errors  are  recorded  for 
each  degree  of  freedom,  and  sorted  based  on  their  relative  perturbation  value  from 
the  nominal.  To  explain  this  concept  graphically,  consider  Figure  3.1. 
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For  each  degree  of  freedom,  the  minimum  error  should  occur  where  the  pre¬ 
dicted  and  true  images  are  precisely  aligned.  Intuitively,  any  deviation  from  this 
point  should  cause  the  error  value  to  increase.  The  data  should  approximately  align 
with  a  parabolic  function,  with  the  errors  increasing  the  further  one  perturbs  from 
the  optimal  value.  This  parabolic  shape  is  approximated  as  a  Gaussian  distribution 
that  is  negated.  The  variance  (cr2)  of  this  Gaussian  distribution  is  the  measurement 
uncertainty  for  that  particular  degree  of  freedom.  The  R  matrix  is  formed  by  placing 
all  six  variances  along  its  diagonal.  The  cross  terms  of  the  R  matrix  are  assumed  to 
be  zero  since  there  is  no  direct  correlation  between  uncertainties. 


Figure  3.1:  Example  error  distribution  about  the  best  image 
match.  The  blue  curve  represents  the  typical  behavior  of  a  single 
degree  of  freedom  when  perturbed  about  the  location  where  the 
true  and  predicted  images  best  match.  A  negated  Guassian 
distribution  is  approximated  based  on  the  shape  of  this  curve, 
whose  variance  is  used  as  the  measurement  uncertainty. 
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3. 3  Summary 

This  chapter  has  shown  how  the  relative  navigation  solution  is  calculated  using 
an  Extended  Kalman  Filter.  The  filter  propagates  the  system  dynamics  over  time 
using  INS  mechanization  for  the  wing  aircraft  and  first  order  Gauss-Markov  approx¬ 
imations  for  the  lead  aircraft.  Using  a  predictive  rendering  approach,  an  estimate  of 
the  lead  aircraft  image  as  seen  from  the  wing  onboard  camera  is  generated  at  the  time 
epochs  where  a  measurement  is  taken.  Through  image  processing  and  perturbation 
techniques,  the  image  update  function  determines  the  relative  position  and  transla¬ 
tion  error  between  the  predicted  and  true  images.  This  serves  as  a  measurement  that 
is  fed  back  to  the  EKF,  which  updates  the  state  estimate.  The  filter  is  then  ready  to 
propagate  forward  again  and  iterate  the  process. 

Based  on  this  development,  Chapter  IV  will  present  analysis  of  the  various  image 
processing  techniques  discussed  in  this  chapter,  as  well  as  multiple  simulations  to  test 
the  relative  navigation  solution  under  various  conditions. 
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IV.  Algorithm  Simulation  and  Analysis 


This  chapter  has  two  major  focuses.  The  first  portion  analyzes  the  three  im¬ 
age  processing  techniques  discussed  earlier  through  numerous  image  processing 
experiments.  The  second  portion  presents  a  performance  analysis  of  the  Extended 
Kalman  Filter  (EKF)  developed  in  Chapter  Iff  using  simulated  trajectories  and  im¬ 
ages. 

4-1  Analysis  of  Image  Processing  Techniques 

This  section  presents  an  analysis  of  the  three  image  processing  techniques  (sum- 
squared  difference,  magnitude  of  gradient,  and  magnitude  of  gradient  with  threshold) 
discussed  in  Chapter  III.  Simulations  were  run  to  determine  the  observability  of  trans¬ 
lation  and  rotation  errors  between  two  given  images  using  each  comparison  technique. 
The  results  of  these  simulations  will  provide  insight  into  how  well  each  technique  will 
detect  errors  between  the  predicted  and  true  images.  All  three  techniques  are  analyzed 
at  the  pre-contact  and  contact  positions  of  the  refueling  formation,  using  infrared  im¬ 
ages  from  the  Northrop  Grumman  test  flights  as  truth  images.  The  plots  shown  in 
this  section  can  be  conceptually  thought  of  as  cost  functions.  The  global  minimum 
should  be  the  optimal  prediction  (i.e.,  the  nominal  estimate  where  the  predicted  image 
exactly  matches  the  truth  image). 

4-1.1  Sum-Squared,  Difference.  Considering  each  degree  of  freedom  for  a 
specific  time  epoch,  the  roll,  pitch,  and  yaw  are  perturbed  from  —20°  to  20°  from 
the  nominal,  and  the  x,  y ,  and  z  body  axes  are  perturbed  from  —40  feet  to  40  feet 
from  the  nominal  position.  These  plots  are  shown  in  Figure  4.1.  The  values  in  the 
plots  are  the  sum-squared  difference  (or  error)  between  a  predicted  image  and  the 
truth  image.  One  can  see  that  the  sum-squared  difference  approach  shows  a  smooth 
rise  in  the  error  as  the  perturbations  deviate  away  from  the  nominal  location  (0  on 
each  plot’s  horizontal  axis).  In  certain  degrees  of  freedom,  the  minimum  error  does 
not  precisely  occur  at  zero.  This  is  clue  to  small  discrepancies  in  the  rendered  model 
dimensions  compared  to  the  images  of  the  true  tanker.  For  example,  differences  in 
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(a)  (b) 

Figure  4.1:  (a)  The  predicted  image  is  independently  per¬ 

turbed  in  roll,  pitch,  and  yaw  about  the  nominal  pre-contact 
position.  The  minimum  error  (using  sum-squared  difference) 
should  occur  where  the  images  are  best  correlated. 

(b)  The  predicted  image  is  independently  perturbed  in  the  x, 
y,  and  z  body  axes  about  the  nominal  pre-contact  position. 

The  minimum  error  (using  sum-squared  difference)  should  occur 
where  the  images  are  best  correlated. 

apparent  engine  size  between  the  rendered  and  true  images  can  be  observed  at  some 
attitudes.  A  sample  surface  plot  is  shown  in  Figure  4.2.  This  particular  plot  shows 
the  relationship  between  both  roll  perturbation  and  body-frame  longitudinal  {x  axis) 
perturbation  about  the  nominal  estimate.  Deviations  from  the  nominal  are  observable 
as  the  surface  smoothly  slopes  to  higher  magnitudes  of  error  the  further  the  predicted 
image  moves  from  the  nominal  image. 

The  next  mode  presented  shows  the  correlation  between  perturbations  about  the 
x  axis  and  perturbations  about  pitch  (Figure  4.3).  There  is  a  definite  observability 
issue  in  this  case,  since  the  very  large  valley  would  make  it  difficult  for  the  image 
processing  system  to  isolate  the  minimum  error.  One  final  mode  at  the  pre-contact 
position  that  shows  observability  issues  is  the  z  axis  translation  coupled  with  pitch, 
which  is  shown  in  Figure  4.4.  Again,  the  undulating  valley  through  this  surface  will 
cause  difficulties  in  determining  the  true  minimum. 
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Figure  4.2:  Example  surface  generated  using  sum-squared  dif¬ 
ference  approach  by  perturbing  both  the  translation  on  the  x 
axis  and  the  amount  of  roll  simultaneously  at  pre-contact  posi¬ 
tion. 
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Figure  4.3:  Poor  observability  example  surface  generated  with 
sum-squared  difference  approach  by  perturbing  both  the  x  axis 
translation  and  amount  of  pitch  simultaneously  at  pre-contact 
position.  The  large  valley  makes  observability  of  the  optimal, 
minimum  location  extremely  difficult. 
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Figure  4.4:  Example  surface  generated  with  sum-squared  dif¬ 
ference  approach  by  perturbing  the  z  axis  translation  and  the 
amount  of  pitch  simultaneously  at  pre-contact  position.  The  val¬ 
ley  that  also  exists  here  makes  the  minimum  difficult  to  locate. 


The  additional  modes  for  the  other  coupled  degrees  of  freedom  demonstrate  fa¬ 
vorable  behavior  similar  to  the  mode  of  Figure  4.2,  showing  good  overall  observability 
of  errors  for  this  technique. 

Now  analyzing  the  contact  position  for  the  sum-squared  difference  approach,  two 
surfaces  are  presented  showing  modes  of  interest.  The  surface  in  Figure  4.5  shows  the 
same  perturbation  as  Figure  4.2,  with  similar,  smooth  behavior.  Figure  4.6  shows  the 
same  perturbation  as  Figure  4.3,  again  showing  poor  observability  since  the  majority 
of  the  surface  is  very  flat  with  no  distinct  minimum. 

With  this  analysis  of  the  sum-squared  difference  technique  as  a  baseline,  the 
next  section  will  consider  the  magnitude  of  gradient  approach  and  analyze  its  potential 
error  detection  capability. 
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Figure  4.5:  Example  surface  generated  using  sum-squared  dif¬ 
ference  approach  by  perturbing  both  the  translation  on  the  x 
axis  and  the  amount  of  roll  simultaneously  at  the  contact  posi¬ 
tion. 


Figure  4.6:  Example  surface  generated  using  sum-squared  dif¬ 
ference  approach  by  perturbing  both  the  translation  on  the  x 
axis  and  the  amount  of  pitch  simultaneously  at  the  contact  po¬ 
sition. 
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4-1.2  Magnitude  of  Gradient.  This  section  will  analyze  the  magnitude  of 
gradient  comparison  technique.  Figure  4.7  shows  similar  plots  as  Figure  4.1.  Here, 
one  notices  that  the  error  increases  more  rapidly  than  the  sum-squared  difference 
technique  as  the  perturbations  vary  from  the  nominal  position.  However,  as  one 
perturbs  farther  and  farther  from  the  nominal,  the  error  begins  to  oscillate  rather 
than  smoothly  increase. 

Figure  4.8  shows  a  similar  surface  as  seen  in  Figure  4.2,  but  compares  magnitude 
of  gradient  between  the  two  images.  Here,  one  sees  that  the  error  increases  more 
rapidly  than  the  sum-squared  difference  method  when  the  image  deviates  from  the 
nominal  position.  There  is  a  distinct,  minimum  value  indicating  the  nominal  location. 


(a)  (b) 

Figure  4.7:  (a)  The  predicted  image  is  independently  per¬ 

turbed  in  roll,  pitch,  and  yaw  about  the  nominal  pre-contact 
position.  The  minimum  error  (using  magnitude  of  gradient) 
should  occur  where  the  images  are  best  correlated. 

(b)  The  predicted  image  is  independently  perturbed  in  the  x, 
y,  and  z  body  axes  about  the  nominal  pre-contact  position. 
The  minimum  error  (using  magnitude  of  gradient)  should  oc¬ 
cur  where  the  images  are  best  correlated. 
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Figure  4.8:  Example  surface  generated  using  magnitude  of 

gradient  approach  by  perturbing  both  the  x  axis  translation  and 
amount  of  roll  simultaneously  at  the  pre-contact  position.  The 
sharp  minimum  indicates  good  observability. 

Additional  surfaces  are  shown  in  Figures  4.9-4.11.  Each  of  these  modes  shows 
numerous  local  minima,  thus  making  the  optimal  position  difficult  to  determine.  With 
multiple  modes  showing  this  type  of  behavior,  the  magnitude  of  gradient  approach 
alone  would  not  be  a  robust  enough  technique  for  precise  error  detection. 
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Figure  4.9:  Example  surface  generated  using  magnitude  of 

gradient  approach  by  perturbing  both  the  y  axis  translation  and 
the  amount  of  yaw  simultaneously  at  the  pre-contact  position. 
There  are  many  local  minima,  indicating  poor  observability. 
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Figure  4.10:  Example  surface  generated  using  magnitude  of 
gradient  by  perturbing  both  y  axis  translation  and  the  amount  of 
roll  simultaneously  at  the  pre-contact  position.  Again,  multiple 
minima  cause  observability  problems. 
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Figure  4.11:  Example  surface  generated  using  magnitude  of 
gradient  approach  by  perturbing  both  z  axis  translation  and  the 
amount  of  pitch  simultaneously  at  the  pre-contact  position.  Os¬ 
cillations  throughout  this  mode  cause  minima  across  the  surface 
which  cause  difficulties  in  image  correlation. 

Now  examining  the  magnitude  of  gradient  technique  while  at  the  contact  posi¬ 
tion,  Figure  4.12  shows  rather  poor  observability  to  the  errors,  but  there  is  a  distinct 
minimum  at  the  nominal  location.  If  the  errors  begin  very  near  the  optimal  position, 
this  approach  could  find  the  correct  minimum.  Figure  4.13  mimics  the  behavior  of 
Figures  4.3  and  4.6,  demonstrating  that  x  axis  translation  error  coupled  with  pitch 
error  proves  to  be  difficult  to  detect. 
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Figure  4.12:  Example  surface  generated  using  magnitude  of 

gradient  by  perturbing  both  x  axis  translation  and  the  amount 
of  roll  simultaneously  at  the  contact  position.  There  is  poor 
observability  across  most  of  the  surface. 
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Figure  4.13:  Example  surface  generated  using  magnitude  of 
gradient  approach  by  perturbing  both  x  axis  perturbation  and 
the  amount  of  pitch  simultaneously  at  the  contact  position.  Sim¬ 
ilar  to  Figures  4.3  and  4.6,  there  is  poor  observability  with  this 
mode,  including  a  valley  in  the  surface. 
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4-1.3  Magnitude  of  Gradient  with  Threshold.  This  section  will  analyze 
the  magnitude  of  gradient  with  threshold  comparison  technique.  Figure  4.14  shows 
the  results  when  a  sample  image  was  perturbed  and  compared  to  the  truth  image 
using  the  gradient  with  threshold  technique.  One  can  see  that  the  plots  do  not  have 
oscillations  that  are  quite  as  severe  as  the  standard  gradient  plots  shown  in  Figure  4.7. 
However,  the  oscillations  are  still  present,  causing  local  minima  at  values  significantly 
different  from  the  nominal  position.  Furthermore,  in  certain  degrees  of  freedom  (such 
as  roll,  pitch,  and  x  axis  perturbation),  the  minimum  value  is  actually  skewed  by  a 
few  degrees/feet. 

An  example  surface  generated  from  the  magnitude  of  gradient  with  threshold 
approach  is  shown  in  Figure  4.15.  This  particular  surface  is  not  quite  as  smooth  as  the 
sum-squared  difference  plot  of  the  same  mode,  but  shows  a  sharper  rise  away  from  the 
nominal  position  as  perturbations  are  introduced.  This  surface  also  does  not  possess 
a  distinct  minimum  at  the  nominal  image  location,  again  causing  a  potential  issue 


(a)  (b) 

Figure  4.14:  (a)  The  predicted  image  is  independently  per¬ 

turbed  in  roll,  pitch,  and  yaw  about  the  nominal  pre-contact 
position.  The  minimum  error  (using  magnitude  of  gradient  with 
threshold)  should  occur  where  the  images  are  best  correlated, 
(b)  The  predicted  image  is  independently  perturbed  in  the  x,  y, 
and  z  body  axes  about  the  nominal  pre-contact  position.  The 
minimum  error  (using  magnitude  of  gradient  with  threshold) 
should  occur  where  the  images  are  best  correlated. 
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Figure  4.15:  Example  surface  generated  using  magnitude  of 
gradient  with  threshold  approach  by  perturbing  both  the  x  axis 
translation  and  the  amount  of  roll  simultaneously  at  the  pre¬ 
contact  position.  The  multiple  minima  near  (0,0)  would  make 
the  optimal  location  difficult  to  detect. 

with  observability.  Furthermore,  examining  the  mode  of  x  axis  perturbation  with 
simultaneous  change  in  pitch,  similar  results  appear  as  those  from  the  sum-squared 
difference  technique.  This  surface  is  shown  in  Figure  4.16. 

Finally,  examining  the  magnitude  of  gradient  with  threshold  while  at  the  contact 
position,  the  first  mode  presented  for  comparison  is  the  x  axis  translation  error  coupled 
with  roll  error.  This  surface  is  shown  in  Figure  4.17.  Very  similar  to  Figure  4.8,  there 
is  a  very  distinct  minimum  that  would  easily  be  detected  when  near  the  optimal  image 
location.  However,  when  the  errors  are  larger,  the  undulations  in  the  outer  portions 
of  the  surface  would  make  error  detection  difficult.  Furthermore,  in  Figure  4.18,  the 
x  axis  error  coupled  with  pitch  error  is  once  again  difficult  to  detect.  Similar  to 
previous  plots,  the  undulations  and  large  valley  could  make  the  distinct  minimum 
difficult  to  detect. 
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Figure  4.16:  Example  surface  generated  using  magnitude  of 

gradient  with  threshold  by  perturbing  the  x  axis  translation  and 
amount  of  pitch  simultaneously  at  the  pre-contact  position.  The 
valley  that  exists  will  hinder  error  detection. 
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Figure  4.17:  Example  surface  generated  using  magnitude  of 
gradient  with  threshold  by  perturbing  the  x  axis  translation  and 
amount  of  roll  simultaneously  at  the  contact  position.  There  is 
a  distinct  minimum  when  near  the  optimal  location,  but  large 
errors  would  fall  on  the  oscillating  outer  portion  of  the  surface. 
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Figure  4.18:  Example  surface  generated  using  magnitude  of 

gradient  with  threshold  by  perturbing  both  the  x  axis  transla¬ 
tion  and  amount  of  pitch  simultaneously  at  the  contact  position. 
Observability  issues  exist  due  to  the  oscillations  and  the  valley. 

4-1-4  Analysis.  Although  the  magnitude  of  gradient  technique  generates  a 
distinct  minimum  in  the  surface  of  Figure  4.8,  its  observability  is  actually  poor  in 
many  of  the  other  modes  shown  in  Figures  4.9-4.11.  Furthermore,  if  the  predicted 
image  shows  large  errors  (i.e. ,  one  or  more  translation  or  rotation  degrees  of  freedom 
are  significantly  different  from  the  values  in  the  true  image),  the  magnitude  of  gradient 
could  struggle  to  identify  the  actual  error.  A  short  example  will  help  illustrate  this 
point. 

Considering  the  two  surfaces  presented  in  Figures  4.2  and  4.8,  assume  that  the 
system  predicts  a  translation  in  the  x  axis  that  has  an  error  of  20  feet.  In  the  case  of 
the  sum-squared  difference  correlation,  if  one  perturbs  the  predicted  value  about  this 
point  (i.e.,  considers  points  on  the  surface  in  the  neighborhood  of  this  prediction)  and 
compares  errors  between  the  new  images  and  the  truth  image,  a  smaller  error  will 
definitely  be  detected  and  the  system  can  adjust  its  prediction  so  that  it  more  closely 
matches  the  true  image.  This  process  can  be  repeated  until  the  prediction  eventually 
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moves  in  the  neighborhood  of  close  correlation  (near  the  point  (0,0)  on  the  graph) 
with  the  true  image  (as  explained  in  the  image  update  discussion  of  Section  3.2.2). 
An  observability  issue  could  exist  very  near  the  nominal  point,  however,  since  this 
region  of  the  surface  is  rather  flat  and  the  error  values  vary  quite  slowly. 

With  the  same  conditions  in  the  magnitude  of  gradient  approach,  small  per¬ 
turbations  about  this  point  on  the  surface  will  likely  yield  very  similar  results  to  the 
predicted  point.  In  other  words,  all  of  the  points  in  the  neighborhood  of  this  predicted 
point  will  yield  similar  errors.  Thus,  the  system  will  have  difficultly  not  only  detecting 
that  the  prediction  is  incorrect,  but  also  determining  which  degree (s)  of  freedom  to 
correct.  Considering  the  magnitude  of  gradient  with  threshold  approach,  the  surface 
shows  similar  behavior  to  the  sum-squared  difference  surface,  although  the  errors  ap¬ 
pear  to  increase  more  rapidly  than  in  Figure  4.2.  Although  this  attribute  would  help 
with  error  detection,  there  is  still  ambiguity  around  the  nominal  point  because  this 
portion  of  the  surface  contains  multiple  minima. 

Based  on  the  above  reasoning,  the  most  robust  approach  appears  to  be  a  com¬ 
bination  of  the  sum-squared  differencing  technique  for  coarse  error  detection,  then 
transitioning  to  the  magnitude  of  gradient  approach  for  fine  tuning.  The  sum-squared 
difference  technique  shows  the  fewest  issues  with  detecting  errors  across  the  entire  sur¬ 
face.  There  are  only  two  modes  (x  axis  translation  error  coupled  with  pitch  error  and 
z  axis  translation  error  coupled  with  pitch  error)  that  show  significant  error  detec¬ 
tion  issues.  This  is  true  for  the  receiver  aircraft  in  both  the  contact  and  pre-contact 
formation  positions.  Based  on  the  surfaces  generated  in  the  three  image  processing 
techniques,  there  is  an  evident  observability  issue  with  x  axis  and  pitch. 

Upon  actually  implementing  the  image  update  in  simulations  (to  be  discussed 
in  detail  in  the  next  section),  the  benefits  of  using  this  two-fold  approach  are  negligi¬ 
ble  in  some  cases.  Through  experimentation,  the  standalone  sum  squared  difference 
approach  performs  almost  identically  to  the  approach  where  sum-squared  difference 
is  used  for  coarse  acquisition  with  a  transition  to  magnitude  of  gradient  for  fine  ad- 
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justment.  However,  there  are  still  potential  observability  issues  with  the  mode  of  x 
body-axis  error  coupled  with  pitch  error,  and  using  the  magnitude  of  gradient  during 
fine  acquisition  did  present  an  advantage.  The  sum-squared  difference  standalone 
approach  demonstrated  some  observability  challenges  in  the  simulations,  which  the 
magnitude  of  gradient  helped  to  mitigate  to  some  extent.  Specific  examples  of  this 
trend  will  be  discussed  in  the  next  section. 

4-2  Filter  Analysis  Using  Generated  Data 

In  order  to  validate  the  filter  design  and  image  update  algorithm,  several  dif¬ 
ferent  scenarios  are  presented  in  this  section.  For  each  scenario,  truth  trajectories 
are  generated  for  both  aircraft.  For  each  wing  trajectory,  simulated  inertial  measure¬ 
ments  are  generated  for  use  in  the  inertial  propagation  portion  of  the  filter.  These 
inertial  measurements  are  corrupted  with  noise  levels  consistent  with  the  grade  of 
the  simulated  inertial  sensor.  There  are  two  grades  of  IMU  modeled  in  the  following 
simulations:  a  navigation  grade  Honeywell  H764-G,  and  a  tactical  grade  Honeywell 
HG-1700.  The  respective  IMU  specifications  are  shown  in  Table  4.1. 

Additionally,  based  on  the  relative  position  and  orientation  between  the  two 
aircraft  along  the  trajectory,  renderings  are  made  of  the  tanker  as  it  should  appear 
from  a  digital  camera  on  the  receiver  aircraft.  Recall  that  Figure  2.5  shows  the  typical 
image  taken  by  the  wing  aircraft  at  the  pre-contact  position,  and  Figure  2.7  shows 
a  corresponding  rendering  at  the  pre-contact  position.  To  provide  the  image  update 
function  with  truth  images,  renderings  are  captured  to  simulate  images  from  the  wing 
onboard  camera.  Each  image  is  also  corrupted  with  a  small  amount  of  translation 
and  orientation  noise  in  order  to  cause  very  slight  perturbations.  Even  in  smooth 
air,  all  aircraft  in  flight  have  small  movements,  so  these  small  perturbations  create 
more  realistic  images.  The  specific  details,  results,  and  analysis  for  each  scenario  are 
presented  in  the  sections  that  follow. 
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Tabic  4.1:  IMU  Performance  Parameters 


Parameter  (Units) 

HG-1700 

H764-G 

Sampling  interval  (ms) 

10.0 

3.906 

Gyro  bias  sigma  (pr) 

1.0 

0.0015 

Gyro  bias  time  constant  (hr) 

1 

1 

Angular  random  walk  (^=) 

0.3 

0.002 

Gyro  scalefactor  sigma  (PPM) 

150 

5 

Accel  bias  sigma  (™) 

0.0098 

2.45e-4 

Accel  bias  time  constant  (hr) 

1 

1 

Velocity  random  walk  (j^=) 

0.57 

0.0143 

Accel  scalefactor  sigma  (PPM) 

300 

100 

4-2.1  Simple  Trajectory  at  Pre-contact  Position.  The  first  simulation  con¬ 
sidered  the  receiver  aircraft  at  the  pre-contact  position,  with  both  aircraft  aligned  in 
straight-and-level  flight.  The  filter  is  initialized  with  the  true  initial  conditions.  The 
initial  state  error  covariance  values  for  this  simulation  are  given  in  Table  4.2.  The 
position,  velocity,  and  attitude  initial  uncertainties  are  the  same  for  both  aircraft. 

The  EKF  is  executed  using  the  generated  INS  measurements  for  propagation 
of  the  wing  aircraft  navigation  state.  At  a  1  Hz  rate,  the  a  priori  relative  navigation 
state  is  used  to  render  the  estimated  position  of  the  tanker,  which  is  then  compared  to 
the  current  truth  image.  Using  the  image  update  function  described  in  Chapter  III, 
the  a  posteriori  estimate  is  computed,  and  the  filter  propagates  to  the  next  image 
update  (1  second  later).  In  each  of  the  following  scenarios,  a  certain  formation  position 
and  IMU  grade  is  considered  and  simulated.  A  30  iteration  Monte  Carlo  analysis  is 
presented  for  each. 

4. 2. 1.1  Navigation  Grade  IMU.  As  a  baseline,  the  relative  errors  and 
covariances  are  shown  for  the  simulation  of  free  inertial  in  Figure  4.19.  In  other 
words,  this  represents  the  typical  errors  when  only  relying  on  the  inertial  sensors  and 
no  measurement  updates.  The  majority  of  the  sample  functions  are  closely  aligned 
with  the  mean. 


Tabic  4.2:  Initial  State  Uncertainties  for  Simulations 


Position  Uncertainty  (North,  East,  and  Down) 

52  (m)2 

Velocity  Uncertainty  (North,  East,  and  Down) 

0.12  (^)2 

Roll  (<fi)  and  Pitch  ( 6 )  Uncertainty 

l2  (deg)2 

Heading  (i(j)  Uncertainty 

52  (deg)'2 

Now  using  the  image  update  function,  Figure  4.20  shows  the  error  in  the  tanker 
position  with  respect  to  the  wing  body  frame  and  the  standard  deviation  (±lcr)  for  the 
ensemble.  Therefore,  this  figure  shows  the  difference  between  the  estimated  and  true 
position  of  the  lead  aircraft  with  respect  to  the  wing  body  frame.  Figure  4.21  displays 
the  position  (estimated  minus  true)  measurement  residuals  during  the  simulations. 
Figure  4.22  displays  the  attitude  (estimated  minus  true)  measurement  residuals  during 
the  simulations. 
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Figure  4.19:  The  mean  relative  position  error  (indicated  by 

the  green  line)  for  the  Monte  Carlo  simulation  using  navigation 
grade  IMU  and  no  image  updates,  expressed  in  the  wing  body 
frame  while  at  the  pre-contact  position.  The  sample  functions 
for  the  30  iterations  are  indicated  by  the  blue,  dotted  lines. 
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Figure  4.20:  Relative  position  error  for  the  Monte  Carlo  sim¬ 
ulation  using  navigation  grade  IMU,  expressed  in  the  wing  body 
frame  while  at  the  pre-contact  position.  The  relative  error  sam¬ 
ple  functions  are  indicated  by  the  blue,  dotted  lines.  The  solid 
green  line  is  the  ensemble  mean,  and  the  red  lines  indicate  the 
ensemble  standard  deviation  (mean  ±lcr  bounds). 

The  plots  indicate  the  filter  can  successfully  track  the  location  of  the  tanker  for 
this  simple  trajectory.  The  mean  y  and  21  errors  in  the  wing  body  frame  remain  under 
one  meter  for  the  duration  of  the  simulation,  while  the  mean  x  errors  remain  under  ap¬ 
proximately  two  meters.  The  measurement  residuals  relative  to  the  navigation  frame 
show  similar  results.  Though  the  sample  functions  of  Figure  4.21  show  more  variance 
in  the  East  and  Down  directions,  the  mean  errors  are  always  less  than  one  meter. 
This  is  consistent  with  the  body  frame  errors  presented  in  Figure  4.20.  The  attitude 
residuals  show  favorable  results,  with  the  majority  of  mean  residual  values  remaining 
less  than  one  degree.  As  demonstrated  in  the  previous  section,  the  image  update  has 
some  observability  issues  in  the  x  axis.  This  is  evident  in  this  simulation  since  the 
x  axis  errors  are  larger  than  the  other  axes. 
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Figure  4.21:  Position  measurement  residuals  for  the  Monte 

Carlo  simulation  using  navigation  grade  IMU  at  the  pre-contact 
position.  The  sample  functions  are  indicated  by  the  blue,  dotted 
lines.  The  ensemble  mean  is  indicated  by  the  green  line. 


- 

0 

10 

20 

30 

40 

50 

60 

0 

10 

20 

30 

40 

50 

60 

/Ns- 

•5 

0  10  20  30  40  50  60 

Elapsed  Time  (s) 


Figure  4.22:  Attitude  measurement  residuals  for  the  Monte 
Carlo  simulation  using  navigation  grade  IMU  at  the  pre-contact 
position.  The  sample  functions  are  indicated  by  the  blue,  dotted 
lines.  The  ensemble  mean  is  indicated  by  the  green  line. 
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4-2. 1.2  Tactical  Grade  IMU.  As  a  performance  comparison,  another 
Monte  Carlo  simulation  of  the  same  trajectory  is  presented  using  the  tactical  grade 
IMU  parameters.  First,  the  free  inertial  simulation  results  are  shown  in  Figure  4.23. 
Secondly,  the  results  for  this  simulation  utilizing  the  image  update  are  shown  in 
Figures  4.24-4.26. 

These  figures  show  similar  performance  to  the  navigation  grade  inertial  sensors 
in  the  y  and  z  axes.  The  errors  in  the  x  body-frame  are  very  consistent  with  those 
in  the  navigation  grade  simulation.  The  measurement  residuals  are  very  similar  in 
behavior  and  magnitude  between  the  tactical  and  navigation  grade  sensors.  For  this 
particular  simulation,  either  grade  of  sensor  yields  accurate  relative  positioning. 
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Figure  4.23:  The  mean  relative  position  error  (indicated  by 

the  green  line)  for  the  Monte  Carlo  simulation  using  tactical 
grade  IMU  and  no  image  updates,  expressed  in  the  wing  body 
frame  while  at  the  pre-contact  position.  The  sample  functions 
for  the  30  iterations  are  indicated  by  the  blue,  dotted  lines. 
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Figure  4.24:  Relative  position  error  for  the  Monte  Carlo  sim¬ 
ulation  using  tactical  grade  IMU,  expressed  in  the  wing  body 
frame  while  at  the  pre-contact  position.  The  relative  error  sam¬ 
ple  functions  are  indicated  by  the  blue,  dotted  lines.  The  solid 
green  curve  is  the  ensemble  mean,  and  the  red  curves  indicate 
the  ensemble  standard  deviation  (mean  ±1  a  bounds). 


4-2.2  Simple  Trajectory  at  Pre- contact  Position  with  Initial  Condition  Errors. 

In  order  to  test  the  filter’s  capability  to  reject  estimation  errors,  another  simulation 
analyzes  significant  errors  in  the  initial  conditions.  In  this  simulation,  the  initial 
position  of  the  lead  aircraft  contains  an  intentional,  random  error  of  5  meters  (lcr)  in 
the  East  direction,  as  well  as  a  roll  error  of  2°  (lcr).  This  is  a  northward  trajectory 
(just  as  the  previous  simulation),  so  the  position  errors  will  appear  in  the  y  axis  of  the 
wing  body  frame.  Due  to  the  intentional  offset  in  East  position  and  roll,  the  initial 
covariance  values  for  these  states  are  increased  to  52  (nr)2  and  22  (deg)2,  respectively. 
The  remaining  initial  conditions  and  covariance  values  are  the  same  as  described  in 
the  previous  simulation  (Table  4.2).  This  30  iteration  Monte  Carlo  simulation  uses  a 
navigation  grade  IMU.  Figures  4.27-4.29  display  the  results  of  this  simulation. 
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Figure  4.25:  Position  measurement  residuals  for  the  Monte 

Carlo  simulation  using  tactical  grade  IMU  at  the  pre-contact 
position.  The  sample  functions  are  indicated  by  the  blue,  dotted 
lines.  The  ensemble  mean  is  indicated  by  the  green  line. 


Figure  4.26:  Attitude  measurement  residuals  for  the  Monte 
Carlo  simulation  using  tactical  grade  IMU  at  the  pre-contact 
position.  The  sample  functions  are  indicated  by  the  blue,  dotted 
lines.  The  ensemble  mean  is  indicated  by  the  green  line. 
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Figure  4.27:  Relative  position  error  for  the  Monte  Carlo  simu¬ 
lation,  expressed  in  the  wing  body  frame  with  intentional  errors 
in  initial  conditions.  The  relative  error  sample  functions  are 
indicated  by  the  bine,  dotted  lines.  The  solid  green  curve  is 
the  ensemble  mean,  and  the  red  curves  indicate  the  ensemble 
standard  deviation  (mean  ±1<t  bounds). 

From  these  figures,  the  filter  demonstrates  its  capability  to  recover  from  initial 
error  estimates.  Within  10  seconds,  the  filter  corrected  for  the  y  body  error,  and 
the  roll  measurement  residuals  were  near  zero.  Again,  the  largest  relative  position 
errors  are  seen  in  the  x  axis  due  to  the  observability  issues.  After  correcting  the 
initial  error,  the  position  residuals  (Figure  4.28)  are  less  than  |1  m|  and  the  attitude 
residuals  (Figure  4.29)  never  exceed  \2°\.  Overall,  the  error  magnitudes  are  similar 
to  the  simulation  that  was  initialized  with  the  true  conditions.  This  is  an  important 
feature  of  the  filter’s  robust  operation,  since  the  true  initial  navigation  state  will 
seldom  be  known  when  executing  a  real  operation. 
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Figure  4.28:  Position  measurement  residuals  for  the  Monte 

Carlo  simulation  with  intentional  errors  in  initial  conditions. 
The  sample  functions  are  indicated  by  the  blue,  dotted  lines. 
The  ensemble  mean  is  indicated  by  the  green  line. 
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Figure  4.29:  Attitude  measurement  residuals  for  the  Monte 
Carlo  simulation  with  intentional  errors  in  initial  conditions. 
The  sample  functions  are  indicated  by  the  blue,  dotted  lines. 
The  ensemble  mean  is  indicated  by  the  green  line. 
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4-2.3  Simple  Trajectory  at  the  Contact  Position.  The  next  simulations 
consider  the  receiver  in  the  contact  position.  The  filter  is  initialized  with  the  true 
initial  conditions,  as  well  as  the  state  error  covariance  values  from  Table  4.2. 

4-2.3. 1  Navigation  Grade  IMU.  The  first  simulation  considers  the 
navigation  grade  IMU.  Figures  4.30-4.32  display  the  results  of  this  simulation.  The 
x  axis  errors  in  the  wing  body  frame  again  show  larger  magnitudes  than  the  y  and  z 
axes.  Recall  from  Figure  2.5  that  only  a  portion  of  the  lead  aircraft  is  visible  at  the 
contact  position.  The  wing  aircraft  camera  field  of  view  can  typically  only  capture 
a  portion  of  the  fuselage  and  the  wings  just  beyond  the  two  inboard  engines.  This 
narrow  view  of  the  tanker  makes  errors  in  x  axis  translation  and  pitch  difficult  to 
detect. 
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Figure  4.30:  Relative  position  error  for  the  Monte  Carlo  simu¬ 
lation  using  navigation  grade  IMU,  expressed  in  the  wing  body 
frame  while  at  the  contact  position.  The  relative  error  sam¬ 
ple  functions  are  indicated  by  the  blue,  dotted  lines.  The  solid 
green  curve  is  the  ensemble  mean,  and  the  red  curves  indicate 
the  ensemble  standard  deviation  (mean  ±lcr  bounds). 
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Figure  4.31:  Position  measurement  residuals  for  the  Monte 

Carlo  simulation  using  navigation  grade  IMU  at  the  contact  po¬ 
sition.  The  sample  functions  are  indicated  by  the  blue,  dotted 
lines.  The  ensemble  mean  is  indicated  by  the  green  line. 
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Figure  4.32:  Attitude  measurement  residuals  for  the  simula¬ 
tion  using  navigation  grade  IMU  at  the  contact  position.  The 
sample  functions  are  indicated  by  the  blue,  dotted  lines.  The 
ensemble  mean  is  indicated  by  the  green  line. 
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With  the  limited  amount  of  the  lead  aircraft  in  the  held  of  view,  these  errors 
would  be  even  more  pronounced  than  in  the  pre-contact  position.  Examining  the 
measurement  residuals,  the  magnitudes  are  very  small  during  the  first  half  of  the 
simulation,  then  begin  to  slightly  increase  in  the  latter  half.  This  is  consistent  with 
the  relative  position  errors  of  Figure  4.30,  where  the  error  and  covariance  in  the  x  body 
axis  suddenly  increase.  The  most  likely  explanation  is  that  a  few  of  the  trajectory 
runs  had  a  bad  measurement  update  at  this  point  in  the  simulation,  consistent  with 
the  decreased  observability  to  the  x  axis  errors. 

4 -2. 3. 2  Tactical  Grade  IMU.  This  same  contact  position  trajectory  is 
simulated  using  a  tactical  grade  IMU.  The  results  are  displayed  in  Figures  4.33-4.35. 
Again,  the  results  here  are  fairly  similar  to  the  previous  tactical  grade  simulation. 
The  tactical  grade  IMU  shows  the  largest  magnitude  of  error  in  the  x  body  axis,  as 
expected.  The  largest  errors  occur  near  the  end  of  the  simulation  (after  50  seconds), 
indicating  that  the  lower  grade  inertial  sensor  does  begin  to  affect  performance  over 
time.  Ultimately,  either  choice  of  sensor  can  yield  a  reasonably  accurate  relative 
position  solution  for  this  simple  trajectory,  though  the  navigation  grade  provides  a 
small  performance  enhancement. 
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Figure  4.33:  Relative  position  error  for  the  Monte  Carlo  sim¬ 
ulation  using  a  tactical  grade  IMU,  expressed  in  the  wing  body 
frame  while  at  the  contact  position.  The  relative  error  sam¬ 
ple  functions  are  indicated  by  the  blue,  dotted  lines.  The  solid 
green  curve  is  the  ensemble  mean,  and  the  red  curves  indicate 
the  ensemble  standard  deviation  (mean  ±lcr  bounds). 

4-3  Summary 

This  chapter  has  highlighted  the  strengths  and  weaknesses  of  three  image  pro¬ 
cessing  techniques.  Both  the  standalone  sum-squared  difference  approach  and  the 
approach  coupling  sum-squared  difference  with  magnitude  of  gradient  give  reasonable 
error  detection  performance  in  the  simulations. 

The  simulated  trajectories  show  the  validity  of  using  predictive  rendering  in  an 
image-aided  EKF  for  relative  positioning.  Though  not  yet  at  a  maturity  level  for 
extremely  precise  (centimeter  level)  positioning,  this  research  has  demonstrated  that 
the  vision-aided  approach  is  a  valuable  tool  for  finding  an  accurate  relative  navigation 
solution.  Furthermore,  this  technique  would  be  a  robust  addition  to  systems  primarily 
relying  on  GPS.  More  on  this  technique’s  future  maturation  and  recommendations 
for  continued  research  are  presented  in  Chapter  V. 
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Figure  4.34:  Position  measurement  residuals  for  the  Monte 

Carlo  simulation  of  straight  and  level  flight  using  a  tactical  grade 
IMU  at  the  contact  position.  The  sample  functions  are  indicated 
by  the  blue,  dotted  lines.  The  ensemble  mean  is  indicated  by 
the  green  line. 


Figure  4.35:  Attitude  measurement  residuals  for  the  Monte 
Carlo  simulation  of  straight  and  level  flight  using  a  tactical  grade 
IMU  at  the  contact  position.  The  sample  functions  are  indicated 
by  the  blue,  dotted  lines.  The  ensemble  mean  is  indicated  by 
the  green  line. 
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V.  Conclusions  and  Recommendations 


Based  on  the  results  from  Chapter  IV,  this  chapter  will  highlight  the  overall 
conclusions  from  this  research  as  well  as  recommendations  for  future  research 
related  to  this  topic. 

5.1  Research  Conclusions 

The  image  processing  analysis  verified  that  the  predictive  rendering  approach 
is  valid  for  detecting  errors  between  a  predicted  and  true  image.  Though  there  are 
several  modes  that  show  observability  issues  (in  particular  the  mode  that  couples 
errors  in  the  x  body  frame  with  errors  in  pitch),  the  simulations  demonstrated  that  the 
image  update  can  still  adequately  detect  residual  errors.  The  approach  of  first  using 
sum-squared  difference  for  coarse  error  detection,  then  transitioning  to  magnitude  of 
gradient  performed  the  best  considering  all  image  processing  simulations. 

Examining  the  simulation  results,  one  can  see  that  the  system  demonstrates 
reasonably  accurate  relative  positioning  during  the  pre-contact  positions.  The  mean 
error  in  relative  position  is  often  less  than  two  meters  with  the  navigation  grade  IMU, 
with  a  few  x  body  frame  errors  exceeding  this  threshold.  The  mean  y  and  z  estimates 
are  nearly  always  within  one  meter  of  the  true  value.  For  the  pre-contact  position, 
this  level  of  accuracy  is  acceptable,  since  there  is  no  physical  contact  between  aircraft 
and  no  danger  of  collision. 

In  the  case  of  initial  condition  errors,  the  filter  rapidly  corrects  these  initial 
errors  and  shows  similar  performance  to  the  previous  simulation.  The  mean  x  axis 
errors  are  the  largest  in  magnitude,  and  generally  stay  less  than  approximately  two 
meters. 

Now  examining  the  contact  position,  recall  from  Chapter  II  the  previous  research 
conducted  by  Ross  (Section  2.10.1).  He  noted  that  an  optical  system  would  need 
to  maintain  accuracy  to  within  five  feet  of  the  nominal  position,  thus  allowing  the 
boom  operator  to  correct  the  residual  error.  Therefore,  an  accuracy  threshold  of 
approximately  1.5  meters  is  appropriate.  The  simulations  for  both  the  navigation 
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and  tactical  grade  inertial  sensors  exceeded  this  threshold  (the  mean  error  deviating 
to  as  much  as  four  meters).  The  current  navigation  solution  alone  would  not  provide 
accurate  enough  positioning  at  the  contact  position  to  yield  satisfactory  and  safe 
operation.  Additional  refinement  is  needed  to  fine  tune  the  accuracy  to  one  meter  or 
better. 

5.2  Recommendations  for  Future  Research 

5.2.1  Additional  Sensor  Analysis.  This  research  used  predictive  rendering 
to  generate  a  measurement  update  during  simple  trajectory  simulations.  Consider¬ 
ing  the  pre-contact  and  contact  positions,  both  navigation  and  tactical  grade  inertial 
sensors  were  simulated  through  various  trajectories.  These  simulations  were  designed 
to  be  completely  passive  for  the  operation  of  the  tanker.  In  other  words,  the  system 
under  test  did  not  rely  on  any  modifications  to  the  tanker  aircraft  (such  as  datalink 
capability  or  optical  markers).  Though  providing  the  baseline  validation  of  the  pre¬ 
dictive  rendering  approach  with  a  completely  passive  system,  there  is  more  research 
to  investigate  with  this  technique. 

Additional  sensor  parameters  (including  ranging  sensors,  datalink  between  air¬ 
craft,  and  dual  GPS  updates)  should  be  analyzed  to  completely  characterize  the  au¬ 
tonomous  refueling  problem.  The  use  of  a  datalink  providing  inertial  or  GPS  position 
updates  would  definitely  be  worth  investigating,  though  this  would  require  tanker  fleet 
modifications.  Furthermore,  including  a  ranging  sensor  on  the  receiver  aircraft  could 
provide  additional  information  to  complement  the  image  update  function  and  make 
it  more  accurate.  An  overall  sensitivity  analysis  is  needed  to  truly  characterize  the 
most  robust,  fault-tolerant  image-aided  system.  Also,  turning  trajectories  are  typical 
in  aerial  refueling  maneuvers  and  should  be  analyzed  using  this  technique. 

5.2.2  Modification  for  Real-Time  Operation.  The  Mat  lab®  simulation  using 
VRML  to  render  the  images  is  a  valid,  yet  time-consuming  approach.  There  is  a 
definite  delay  between  commanding  the  rendering  parameters  and  the  VR  window 
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actually  displaying  the  rendering,  even  noticeable  to  the  human  eye.  In  order  to 
mitigate  this,  other  rendering  methods  should  be  explored.  Specifically,  directly  using 
the  graphics  processing  unit  to  render  would  likely  be  the  optimal  approach  to  achieve 
high  enough  frame  rates  to  render  and  process  many  images  every  second. 

This  research  used  an  EKF,  which  generates  only  one  navigation  state  estimate. 
Thus,  perturbations  about  this  state  estimate  were  needed  to  determine  the  residual 
errors  between  the  predicted  and  measured  relative  position  and  orientation.  Ren¬ 
dering  each  of  these  perturbations  and  iteratively  adjusting  the  prediction  to  find  the 
closest  matching  image  was  a  lengthy,  time-consuming  process.  For  real-time  opera¬ 
tion,  this  approach  would  not  be  practical.  To  enhance  the  speed  of  the  image  update, 
an  Unscented  Kalman  Filter  (UKF)  should  be  implemented.  With  a  UKF,  a  large 
number  of  sigma  points  would  be  generated  instead  of  a  single  state  estimate.  With 
each  sigma  point  as  a  prediction  (and  thus  predicted  image),  the  image  update  could 
develop  some  type  of  hypersurface  to  characterize  how  the  errors  vary  with  the  sigma 
point  spread.  Based  on  this  distribution  of  error  values,  the  update  function  could 
determine  where  the  minimum  error,  and  thus  optimal  state  estimate,  occurs.  The 
function  would  still  be  generating  a  large  number  of  rendered  images,  but  would  not 
be  performing  the  time-consuming  perturbation  process.  This  would  allow  the  image 
update  to  capture  a  truth  image  from  the  digital  camera,  compare  the  sigma  point 
predictions  against  the  truth,  determine  the  nominal  estimate,  and  pass  the  nominal 
back  to  the  filter  for  propagation  to  the  next  measurement.  All  of  this  must  occur 
quickly  enough  to  capture  truth  images  from  the  digital  camera  at  a  1  Hz  or  possibly 
2  Hz  rate. 

Finally,  with  a  filter  and  image  update  configuration  that  can  operate  in  real¬ 
time,  a  controller  is  needed  to  correct  the  receiver  aircraft  position  based  on  the 
navigation  state  estimate.  This  could  be  a  simple  regulator  that  drives  the  error 
state  to  zero,  thus  keeping  the  receiver  aircraft  at  the  nominal  pre-contact  or  contact 
position,  as  appropriate. 
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